#include <stdlib.h>
#include "LegsMakerPiPlusAsymJul00.h"
#include "LegsRun.h"

ClassImp(LegsMakerPiPlusAsymJul00) 
//____________________________________________________________________

 LegsMakerPiPlusAsymJul00::LegsMakerPiPlusAsymJul00(Int_t add)
    : LegsMaker ("PiPlusAsymJul00","Pi+ Asymmetry in XT",  add)
{
  // constructor of the class.  Data member "it" initialized to be
  // ponter to the LegsRun object which is processed right now. The maker
  // added to the list of makers in LegsRun object, so that it knows
  // about it. subdirectory in output file created for this maker
  // (the directory has the same name as the maker does).
  
  Init(); this->cd("../");
}  

 void LegsMakerPiPlusAsymJul00::Init()
{
  // user initialization, create histos  etc
  char histName[100];
  char *event[5] = {"yld","yld*pol","true","acc","rest"};

  MaxEgy      = 320.0;  // MLUV data (hd)
  TagBins     = 16;
  TrueTagPeak = 35.0;  // MLUV data (hd)
  TagPeakDist = 18.9;
  LowTrueTag  = TrueTagPeak - TagPeakDist / 2.0;
  HighTrueTag = TrueTagPeak + TagPeakDist / 2.0;
  HighAccTag  = TrueTagPeak + 3.5 * TagPeakDist;

  //  accscal = 1./3.;
  accscal = 1./2.76;
  Lambda  = 2./TMath::Pi();

  v1 = new TLorentzVector();

  sprintf(histName,"Tagger_TDC");
  fHistTagTDC = new TH2F(histName,histName,121,0,120,66,0,65);
  // cout << " creating histo " << histName << endl;

  sprintf(histName,"pol");
  fPolHist = new TH1F(histName,histName,36,-1.,8.);
  //  cout << " creating histo  " << histName << endl;
  
  for(UInt_t face=0; face < 4; face++)
    {  
      sprintf(histName,"XV_ADC_face%i",face+1);
      fHistXvAdc[face] = new TH1F(histName,histName,100,0,256);
      // cout << " creating histo " << histName << endl;
    } 

  for(Int_t k=0; k<5; k++)
    {
      sprintf(histName,"EdE_%s",event[k]);
      fHistPiPlusAsymEdE[k][0][0] = new TH2F(histName,histName,100,0,350,100,0,10);
      //  cout << " creating histo  " << histName << endl;

      for(Int_t pol=0; pol < 3; pol++) // number of pol states
        {
          sprintf(histName,"Eg_pol%i_%s",pol,event[k]);
          fHistPiPlusAsym1Eg[k][pol] = new TH1F(histName,histName,70,0.0,420.0);
	  fHistPiPlusAsym1Eg[k][pol] ->Sumw2();
          //	cout << " creating histo  " << histName << endl;
          sprintf(histName,"Theta_pol%i_%s",pol,event[k]);
          fHistPiPlusAsym1Theta[k][pol] = new TH1F(histName,histName,18,0.0,180.0);
	  fHistPiPlusAsym1Theta[k][pol] ->Sumw2();
          //	cout << " creating histo  " << histName << endl;
          sprintf(histName,"Phi_pol%i_%s",pol,event[k]);
          fHistPiPlusAsym1Phi[k][pol] = new TH1F(histName,histName,37,-5.0,365.0);
	  fHistPiPlusAsym1Phi[k][pol] ->Sumw2();
          //	cout << " creating histo  " << histName << endl;
	  sprintf(histName,"Eg_vs_phi_pol%i_%s",pol,event[k]);
	  fHistPiPlusAsym2Egphi[k][pol] = new TH2F(histName,histName,32,228.0,420.0,37,-5.0,365.0);
	  fHistPiPlusAsym2Egphi[k][pol] ->Sumw2();
	  //  cout << " creating histo  " << histName << endl;
	  sprintf(histName,"Eg_vs_theta_pol%i_%s",pol,event[k]);
	  fHistPiPlusAsym2Egtheta[k][pol] = new TH2F(histName,histName,32,228.0,420.0,18,0.0,180.0);
	  fHistPiPlusAsym2Egtheta[k][pol] ->Sumw2();
	  //  cout << " creating histo  " << histName << endl;
	  sprintf(histName,"Theta_vs_phi_pol%i_%s",pol,event[k]);
	  fHistPiPlusAsym2Thetaphi[k][pol] = new TH2F(histName,histName,18,0.0,180.0,37,-5.0,365.0);
	  fHistPiPlusAsym2Thetaphi[k][pol] ->Sumw2();
	  //  cout << " creating histo  " << histName << endl;
	  sprintf(histName,"HR_Theta_vs_phi_pol%i_%s",pol,event[k]);
	  fHistPiPlusAsymHRThetaphi[k][pol] = new TH2F(histName,histName,180,0.0,180.0,360,0.0,360.0);
	  fHistPiPlusAsymHRThetaphi[k][pol] ->Sumw2();
	  //  cout << " creating histo  " << histName << endl;
	  sprintf(histName,"HR_CosTheta_vs_phi_pol%i_%s",pol,event[k]);
	  fHistPiPlusAsymHRCosThetaphi[k][pol] = new TH2F(histName,histName,100,-1.0,1.0,360,0.0,360.0);
	  fHistPiPlusAsymHRCosThetaphi[k][pol] ->Sumw2();
	  //  cout << " creating histo  " << histName << endl;

          for(Int_t face=0; face < 4; face++)
            {
              sprintf(histName,"Eg_pol%i_face%i_%s",pol,face+1,event[k]);
              fHistPiPlusAsymEgface[k][pol][face] = new TH1F(histName,histName,70,0.0,420.0);
	      fHistPiPlusAsymEgface[k][pol][face] ->Sumw2();
	      //  cout << " creating histo  " << histName << endl;
              sprintf(histName,"Theta_pol%i_face%i_%s",pol,face+1,event[k]);
              fHistPiPlusAsymThetaface[k][pol][face] = new TH1F(histName,histName,18,0.0,180.0);
	      fHistPiPlusAsymThetaface[k][pol][face] ->Sumw2();
	      //  cout << " creating histo  " << histName << endl;
	      sprintf(histName,"Eg_vs_theta_pol%i_face%i_%s",pol,face+1,event[k]);
	      fHistPiPlusAsym2EgthetaLBface[k][pol][face] = new TH2F(histName,histName,16,228.0,420.0,18,0.0,180.0);
	      fHistPiPlusAsym2EgthetaLBface[k][pol][face] ->Sumw2();
	      //  cout << " creating histo  " << histName << endl;
	      sprintf(histName,"Eg_vs_thetacm_pol%i_face%i_%s",pol,face+1,event[k]);
	      fHistPiPlusAsym2EgthetaCMface[k][pol][face] = new TH2F(histName,histName,16,228.0,420.0,18,0.0,180.0);
	      fHistPiPlusAsym2EgthetaCMface[k][pol][face] ->Sumw2();
	      //  cout << " creating histo  " << histName << endl;
              for(Int_t tagChannel=0; tagChannel< TagBins; tagChannel++)
                {
	          sprintf(histName,"ThetaCM_pol%i_face%i_tag%i_%s",pol,face+1,tagChannel+1,event[k]);
	          fHistPiPlusAsym1ThetaCMface[k][pol][face][tagChannel] = new TH1F(histName,histName,18,0.0,180.0);
	          fHistPiPlusAsym1ThetaCMface[k][pol][face][tagChannel] ->Sumw2();
	      //  cout << " creating histo  " << histName << endl;
	        }
	    }
        }

      for(Int_t face=0; face < 4; face++)
        {
          sprintf(histName,"EdE%i_%s",face+1,event[k]);
          fHistPiPlusAsymEdEdiffFaces[k][face] = new TH2F(histName,histName,100,0,350,100,0,10); 
          //  cout << " creating histo  " << histName << endl;
        }
  
      for(Int_t tagChannel=0; tagChannel< TagBins; tagChannel++)
        {
          sprintf(histName,"EdE_tag%i_%s",tagChannel+1,event[k]);
          fHistPiPlusAsymEdEdiffTagChannel[k][tagChannel] = new TH2F(histName,histName,100,0,350,100,0,10);
          //      cout << " creating histo  " << histName << endl;
	  sprintf(histName,"Epion_vs_Phi_tag%i_%s",tagChannel+1,event[k]);
	  fHistPiPlusAsym2PiPhi[k][tagChannel] = new TH2F(histName,histName,60,0,300,37,-5.0,365.0);
	  //	  cout << " creating histo  " << histName << endl;
	  sprintf(histName,"Epion_vs_ThePhi_tag%i_%s",tagChannel+1,event[k]);
	  fHistPiPlusAsym3PiThePhi[k][tagChannel] = new TH3F(histName,histName,50,0,250,18,0.0,180.0,37,-5.0,365.0);
	  //	  cout << " creating histo  " << histName << endl;

          for(Int_t face=0; face < 4; face++)
            {
	      sprintf(histName,"EdE%i_tag%i_%s",face+1,tagChannel+1,event[k]);
	      fHistPiPlusAsym1[k][face][tagChannel] = new TH2F(histName,histName,100,0,350,100,0,10);
	      //	cout << " creating histo  " << histName << endl;
	      sprintf(histName,"Epion%i_tag%i_%s",face+1,tagChannel+1,event[k]);
	      fHistPiPlusAsym1Pion[k][face][tagChannel] = new TH1F(histName,histName,351,0,350);
	      //	cout << " creating histo  " << histName << endl;
	      sprintf(histName,"Epion_vs_Theta_face%i_tag%i_%s",face+1,tagChannel+1,event[k]);
	      fHistPiPlusAsym2PiTheta[k][face][tagChannel] = new TH2F(histName,histName,60,0,300,18,0.0,180.0);
	      //	cout << " creating histo  " << histName << endl;
            }

          for(Int_t pol=0; pol < 3; pol++) // number of pol states
            {
	      sprintf(histName,"angles_pol%i_tag%i_%s",pol,tagChannel+1,event[k]);
	      //  cout << " creating histo  " << histName << endl;
    	      fHistPiPlusAsym[k][pol][tagChannel] = new TH2F(histName,histName,37,-5.0,365.0,18,0.0,180.0);
	      fHistPiPlusAsym[k][pol][tagChannel] ->Sumw2();
	      //  fHistPiPlusAsym[k][i][j] = new TH2F(histName,histName,144,-0.1,TWOPI+0.1,40,-0.1,PI+0.1);
            }
	}
    }
}

 bool LegsMakerPiPlusAsymJul00::Make()
{
  // returns TRUE if charged particle in XT and no charged particle in
  // bar vetos.
  // returns FALSE otherwise
  // Makes EdE histos for cutting out pi+.
  
//   cout << " PiPlusAsym 1" << endl;
  Float_t DEtmp;
  bool flag=false; 
  //  char histName[50];
  Int_t pnNpart=0;
  Int_t face;
  Int_t xtFace[4],xvFace[4];
  //  Float_t allFaces=0.0;
  Float_t xvAdc[4];
  Double_t angleSin;
  //  Float_t perpendicular = fXt[1]->GetX(); // for all faces of xt distance to target from center of a face
  // count number of hits per face in XT and XV
  UInt_t i;
  UInt_t hid;
  for(hid=0; hid<4; hid++)
    {
      xtFace[hid] = 0;
      xvFace[hid] = 0;
      xvAdc[hid]  = 0;
    }
  
  // make sure there is something in XT and XV and nothing in PV
  
  if(it->rxt_count <=0 || it->rpv_count > 0 || it->rxv_count == 0) return false;
  //  if(it->rxt_count <=0 || it->rpv_count > 0 || it->rxv_count == 0  || it->GetTagCount() != 1) return false;
  //  if(it->rxt_count <=0 || it->rpv_count > 0 || it->rxv_count == 0  || it->cpn_count == 0 || it->GetTagCount() != 1) return false;
  
  if(it->GetTagID(0) <=4) return false; // bad channels
  Int_t SumBins = 60 / (TagBins-1);
  Int_t tagChannel=(it->GetTagID(0)-SumBins-1) / SumBins + 1;
  //  Int_t tagChannel=(it->GetTagID(0)-5)/5+1;
  //  cout << it->GetTagID(0) << " " << tagChannel<< endl;
  
  // fill in arrays with number of hits in a face and ADC for a face for XV
  
  Int_t itr;
  
  for(itr = 0; itr < it->GetCount();itr++)
    {
      if(it->GetDet(itr) == 1) // XV
	{
	  xvFace[(it->GetFace(itr))-1]++;
	  xvAdc[(it->GetFace(itr))-1] = it->GetAdc1(itr);
	  //	  allFaces=allFaces+1;
	}
    }
  
  // loop over all clusters-particles
  Int_t tmpDet;
  for(i=0;i<it->fPartVec.size();i++)
    {
      tmpDet = it->fPartVec[i].GetDet();
      if(tmpDet == 2) // XT
	{
	  face = it->fPartVec[i].GetFace();
	  xtFace[face-1]++;
	}
      else if(tmpDet == 5) // PN
	{
	  pnNpart++;
	}
    }

    static Int_t polstate;
    static Float_t norm1;
    static Float_t yieldpol1;
    static Float_t norm1acc;
    static Float_t yieldpol1acc;

    if (it->GetRegPol() == 0)
      polstate = 0;
    else if (it->GetRegPol() == 3)
      polstate = 1;
    else if (it->GetRegPol() == 4)
      polstate = 2;
    else
      //      return flag;
      polstate = 100;

    assert(polstate >= 0);
    assert(polstate <= 2);
    //    norm  = 1./(it->fDB->GetCg3(0,polstate,"int"));
    norm1 = 1./(it->fDB->GetCg3(it->GetTagID(0),polstate,"int"));
    //    normacc = norm * accscal;
    norm1acc = norm1 * accscal;
    if (polstate == 2)
      {
	//        yieldpol  = Lambda * norm * (it->fDB->GetPhotPol(it->GetTagID(0),2));
        yieldpol1 = Lambda * norm1 * (it->fDB->GetPhotPol(it->GetTagID(0),3));
      }
    else if (polstate == 1)
      {
	//        yieldpol  = Lambda * norm * (it->fDB->GetPhotPol(it->GetTagID(0),1));
        yieldpol1 = Lambda * norm1 * (it->fDB->GetPhotPol(it->GetTagID(0),4));
      }
    else
      {
	//        yieldpol  = Lambda * norm * (it->fDB->GetPhotPol(it->GetTagID(0),0));
        yieldpol1 = Lambda * norm1 * (it->fDB->GetPhotPol(it->GetTagID(0),0));
      }

    //    yieldpolacc = yieldpol * accscal;
    yieldpol1acc = yieldpol1 * accscal;

    //    cout << "I am after polstate " << polstate << "   " << norm << "   " << endl;

//   cout << " PiPlusAsym 2" << endl;
  for(i=0; i < it->fPartVec.size();i++)
    {
      if(it->fPartVec[i].GetDet() == 2)
	{ 
	  // fill histo with EdE
	  face = it->fPartVec[i].GetFace();
	  //	  cout << xtFace[face] << " " << xvFace[face] <<  " face = " << face << endl;
	  //	  if(xtFace[face]==1 && xvFace[face]==1 && ((face == 1 && xvFace[3]==0) || (face == 3 && xvFace[1]==0) || (face == 2 && xvFace[4]==0) || (face == 4 && xvFace[2]==0))) // one cluster per XT face
	  if(xtFace[face-1]==1 && xvFace[face-1]==1 && (face == 1 || face == 3 || face == 2 || face == 4 )) // one cluster per XT face
	    {
	      // make sure it is not cosmics - they  go through and leave signal in opposite face of top and bottom faces
	      
	      // correct for angle dependence of thickness of veto paddle which is passed through by particle
	      //      cout << xtFace[face] << " " << xvFace[face] << endl;
	      flag=true;
	      angleSin = fabs(sin(it->fPartVec[i].GetTheta()));
	      assert(angleSin <= 1);
	      //	      DEtmp = xvAdc[face]*angleSin;
	      DEtmp = it->fDB->GetXvEgy(it->fPartVec[i].GetTheta(),it->fPartVec[i].GetPhi(),(Int_t)xvAdc[face-1],face);	      
	      DEtmp = DEtmp*angleSin;

	      bool xc_only = true;
	      for(UInt_t kf=0; kf < it->fPartVec[i].Nhits(); kf++)
              if(it->fPartVec[i].GetDet() != 3) 
		{
		  xc_only=false;   
		  break;
	        }
	      //if(!xc_only)
	      //{
	      //	      it->Show(it->GetCurrent());
	      Float_t partEgy = it->fPartVec[i].GetLgt();
	      UShort_t tagid  = it->GetTagID(0);
	      Float_t tagtdc = it->fDB->GetTagTof(tagid,it->GetTagTdc(0));
	      //	      cout << it->fPartVec[i].GetLgt() << " " << DEtmp << endl;
	      //	      AddToEventList();
	      
	      if(partEgy + 50.0*DEtmp > 50)
		{  
		  // for 4 faces separately
		  if(it->fPartVec[i].GetTheta() > (kPI/2.0 - 0.1) && it->fPartVec[i].GetTheta() < (kPI/2.0 + 0.1) )
		    {
		      // reverse numbering of faces and phi ...
		      bool phiYes = false;
		      Float_t tPhi = it->fPartVec[i].GetPhi();
		      switch (face)
			{
			case 1:
			  if(tPhi > (kTWOPI-0.1)  || tPhi < 0.1)
			    phiYes = true;
			  break;
			case 2:
			  if(tPhi > (3.0*kPI/2.0-0.1)  &&  tPhi < (3.0*kPI/2.0 + 0.1))
			    phiYes = true;
			  break;
			case 3:
			  if(tPhi > (kPI-0.1)  &&  tPhi < (kPI + 0.1))
			    phiYes = true;
			  break;
			case 4:
			  if(tPhi > (kPI/2.0-0.1)  &&  tPhi < (kPI/2.0 + 0.1))
			    phiYes = true;
			  break;
			default:
			  assert("face number must be in range 1 to 4"=="!?");
			}
		      if(phiYes)
			{
			  assert(face>=1);
			  assert(face<=4);
   			  fHistXvAdc[face-1]->Fill(xvAdc[face-1]);
			}
		    }
		  
		  // Tagger TDC vs. Tagger ID
   		  fHistTagTDC->Fill(tagtdc,tagid);

		  assert(face >= 1);
		  assert(face <= 4);
		  assert(tagChannel >= 0);
		  assert(tagChannel <  TagBins);
		  // for 4 faces and  12 tag chanels separately
		  if(partEgy + 50.0*DEtmp > 50)
		    {  
		      if (tagtdc >= LowTrueTag && tagtdc < HighTrueTag)
			{
		          fHistPiPlusAsym1[0][face-1][tagChannel]->Fill(partEgy,DEtmp,1);
		          fHistPiPlusAsym1[1][face-1][tagChannel]->Fill(partEgy,DEtmp,1);
		          fHistPiPlusAsym1[2][face-1][tagChannel]->Fill(partEgy,DEtmp);
			}
		      else if (tagtdc >= HighTrueTag && tagtdc < HighAccTag)
			{
		          fHistPiPlusAsym1[0][face-1][tagChannel]->Fill(partEgy,DEtmp,-accscal);
		          fHistPiPlusAsym1[1][face-1][tagChannel]->Fill(partEgy,DEtmp,-accscal);
		          fHistPiPlusAsym1[3][face-1][tagChannel]->Fill(partEgy,DEtmp);
			}
		      else
		        fHistPiPlusAsym1[4][face-1][tagChannel]->Fill(partEgy,DEtmp);
		    }
		  // for 4 faces separately and all tag chanels together
		  if( (fHistPiPlusAsymEdEdiffFaces[0][face-1] != NULL) ||
		      (fHistPiPlusAsymEdEdiffFaces[1][face-1] != NULL) ||
		      (fHistPiPlusAsymEdEdiffFaces[2][face-1] != NULL) ||
		      (fHistPiPlusAsymEdEdiffFaces[3][face-1] != NULL) ||
		      (fHistPiPlusAsymEdEdiffFaces[4][face-1] != NULL) )
		    {  
		      if (tagtdc >= LowTrueTag && tagtdc < HighTrueTag)
			{
		          fHistPiPlusAsymEdEdiffFaces[0][face-1]->Fill(partEgy,DEtmp,1);
		          fHistPiPlusAsymEdEdiffFaces[1][face-1]->Fill(partEgy,DEtmp,1);
		          fHistPiPlusAsymEdEdiffFaces[2][face-1]->Fill(partEgy,DEtmp);
			}
		      else if (tagtdc >= HighTrueTag && tagtdc < HighAccTag)
			{
		          fHistPiPlusAsymEdEdiffFaces[0][face-1]->Fill(partEgy,DEtmp,-accscal);
		          fHistPiPlusAsymEdEdiffFaces[1][face-1]->Fill(partEgy,DEtmp,-accscal);
		          fHistPiPlusAsymEdEdiffFaces[3][face-1]->Fill(partEgy,DEtmp);
			}
		      else
		        fHistPiPlusAsymEdEdiffFaces[4][face-1]->Fill(partEgy,DEtmp);
		    }
		  // SUM all faces but separate tag chanels
		  if( (fHistPiPlusAsymEdEdiffTagChannel[0][tagChannel] != NULL) ||
		      (fHistPiPlusAsymEdEdiffTagChannel[1][tagChannel] != NULL) ||
		      (fHistPiPlusAsymEdEdiffTagChannel[2][tagChannel] != NULL) ||
		      (fHistPiPlusAsymEdEdiffTagChannel[3][tagChannel] != NULL) ||
		      (fHistPiPlusAsymEdEdiffTagChannel[4][tagChannel] != NULL) )
		    {  
		      if (tagtdc >= LowTrueTag && tagtdc < HighTrueTag)
			{
		          fHistPiPlusAsymEdEdiffTagChannel[0][tagChannel]->Fill(partEgy,DEtmp,1);
		          fHistPiPlusAsymEdEdiffTagChannel[1][tagChannel]->Fill(partEgy,DEtmp,1);
		          fHistPiPlusAsymEdEdiffTagChannel[2][tagChannel]->Fill(partEgy,DEtmp);
			}
		      else if (tagtdc >= HighTrueTag && tagtdc < HighAccTag)
			{
		          fHistPiPlusAsymEdEdiffTagChannel[0][tagChannel]->Fill(partEgy,DEtmp,-accscal);
		          fHistPiPlusAsymEdEdiffTagChannel[1][tagChannel]->Fill(partEgy,DEtmp,-accscal);
		          fHistPiPlusAsymEdEdiffTagChannel[3][tagChannel]->Fill(partEgy,DEtmp);
			}
		      else
		        fHistPiPlusAsymEdEdiffTagChannel[4][tagChannel]->Fill(partEgy,DEtmp);
		    }
		  
		  if( (fHistPiPlusAsymEdE[0][0][0] != NULL) ||
		      (fHistPiPlusAsymEdE[1][0][0] != NULL) ||
		      (fHistPiPlusAsymEdE[2][0][0] != NULL) ||
		      (fHistPiPlusAsymEdE[3][0][0] != NULL) ||
		      (fHistPiPlusAsymEdE[4][0][0] != NULL) )
		    {  
		      if (tagtdc >= LowTrueTag && tagtdc < HighTrueTag)
			{
		          fHistPiPlusAsymEdE[0][0][0]->Fill(partEgy,DEtmp,1);
		          fHistPiPlusAsymEdE[1][0][0]->Fill(partEgy,DEtmp,1);
		          fHistPiPlusAsymEdE[2][0][0]->Fill(partEgy,DEtmp);
			}
		      else if (tagtdc >= HighTrueTag && tagtdc < HighAccTag)
			{
		          fHistPiPlusAsymEdE[0][0][0]->Fill(partEgy,DEtmp,-accscal);
		          fHistPiPlusAsymEdE[1][0][0]->Fill(partEgy,DEtmp,-accscal);
		          fHistPiPlusAsymEdE[3][0][0]->Fill(partEgy,DEtmp);
			}
		      else
		        fHistPiPlusAsymEdE[4][0][0]->Fill(partEgy,DEtmp);
		    }
		}
	      //	      if( DEtmp > 0.8 && partEgy+DEtmp*2.0 > 30  && partEgy+80.0*DEtmp < 400) 
	      if( DEtmp > 0.8 && partEgy+DEtmp*2.0 > 20  && partEgy+80.0*DEtmp < 400) 
		{
		  //   cout << fPartVec[i].Nhits() << endl;
		  static Float_t tmpFloat;
		  tmpFloat = it->GetRegPol();
		  fPolHist->Fill(tmpFloat);
		  static Float_t RadPartPhi;
		  static Float_t RadPartTheta;
		  static Float_t PartPhi;
		  static Float_t PartTheta;
		  static Float_t PhotEgy;
		  RadPartPhi   = it->fPartVec[i].GetPhi();
		  RadPartTheta = it->fPartVec[i].GetTheta();
		  PartPhi      = RadPartPhi / TMath::Pi() * 180.;
		  PartTheta    = RadPartTheta / TMath::Pi() * 180.;
		  PhotEgy      = it->fDB->GetTagEgy(tagid);

		  // do a Lorentz-Transformation of the pion from Lab to CM
		  Double_t px,py,pz,E;
		  const Double_t mProt = 938.27231;
		  const Double_t mPion = 139.56995;
		  px = 0.; py = 0.; pz = PhotEgy; E = pz + mProt;
		  v1->SetPxPyPzE(px,py,pz,E);
		  Double_t beta  = v1->Beta();
		  Double_t gamma = v1->Gamma();

		  Double_t PionPxlb,PionPylb,PionPzlb,PionElb,PionPlb,PionTlb,PionPplb,PionPtlb;
		  PionTlb  = partEgy + DEtmp;
		  PionPlb  = sqrt(PionTlb * (PionTlb + 2.0*mPion));
		  PionElb  = sqrt(PionPlb*PionPlb + mPion*mPion);
		  PionPxlb = PionPlb * sin(RadPartTheta) * cos(RadPartPhi);
		  PionPylb = PionPlb * sin(RadPartTheta) * sin(RadPartPhi);
		  PionPzlb = PionPlb * cos(RadPartTheta);
		  PionPplb = PionPzlb;
		  PionPtlb = sqrt(PionPxlb*PionPxlb + PionPylb*PionPylb);

		  Double_t PionPpcm,PionPtcm,PartThetaCM;
		  PionPpcm = gamma * (PionPplb - beta*PionElb);
		  PionPtcm = PionPtlb;
		  PartThetaCM = atan2(PionPtcm,PionPpcm) / TMath::Pi() * 180.;

		  assert(polstate >= 0);
		  assert(polstate <= 2);
		  assert(tagChannel >= 0);
		  assert(tagChannel <  TagBins);

		  //		  if ( !((PartPhi > 70. && PartPhi < 110.) || (PartPhi > 160. && PartPhi < 200.) || 
		  //		         (PartPhi > 250. && PartPhi < 290.) || (PartPhi > 340. && PartPhi <= 360.) ||
		  //		         (PartPhi >= 0. && PartPhi < 20.)) ) return false;
		      if (tagtdc >= LowTrueTag && tagtdc < HighTrueTag)
			{
			   fHistPiPlusAsym1Eg[0][polstate] ->Fill(PhotEgy,norm1);
			   fHistPiPlusAsymEgface[0][polstate][face-1] ->Fill(PhotEgy,norm1);
		           fHistPiPlusAsym2EgthetaLBface[0][polstate][face-1] ->Fill(PhotEgy,PartTheta,norm1);
		           fHistPiPlusAsym2EgthetaCMface[0][polstate][face-1] ->Fill(PhotEgy,PartThetaCM,norm1);
		           fHistPiPlusAsym2Egphi[0][polstate] ->Fill(PhotEgy,PartPhi,norm1);
		           fHistPiPlusAsym2Egtheta[0][polstate] ->Fill(PhotEgy,PartTheta,norm1);
			   if (PhotEgy < MaxEgy)
			     {
		               fHistPiPlusAsym[0][polstate][tagChannel] ->Fill(PartPhi,PartTheta,norm1);
			       fHistPiPlusAsymThetaface[0][polstate][face-1] ->Fill(PartTheta,norm1);
			       fHistPiPlusAsym1Theta[0][polstate] ->Fill(PartTheta,norm1);
			       fHistPiPlusAsym1Phi[0][polstate] ->Fill(PartPhi,norm1);
		               fHistPiPlusAsym2Thetaphi[0][polstate] ->Fill(PartTheta,PartPhi,norm1);
		               fHistPiPlusAsymHRThetaphi[0][polstate] ->Fill(PartTheta,PartPhi,norm1);
		               fHistPiPlusAsymHRCosThetaphi[0][polstate] ->Fill(cos(RadPartTheta),PartPhi,norm1);
			       fHistPiPlusAsym1ThetaCMface[0][polstate][face-1][tagChannel] ->Fill(PartThetaCM,norm1);
			       fHistPiPlusAsym1Pion[0][face-1][tagChannel] ->Fill(partEgy,1);
			       fHistPiPlusAsym2PiTheta[0][face-1][tagChannel] ->Fill(partEgy,PartTheta,1);
			       fHistPiPlusAsym2PiPhi[0][tagChannel] ->Fill(partEgy,PartPhi,1);
			       fHistPiPlusAsym3PiThePhi[0][tagChannel] ->Fill(partEgy,PartTheta,PartPhi,1);
			     }

			   fHistPiPlusAsym1Eg[1][polstate] ->Fill(PhotEgy,yieldpol1);
			   fHistPiPlusAsymEgface[1][polstate][face-1] ->Fill(PhotEgy,yieldpol1);
		           fHistPiPlusAsym2EgthetaLBface[1][polstate][face-1] ->Fill(PhotEgy,PartTheta,yieldpol1);
		           fHistPiPlusAsym2EgthetaCMface[1][polstate][face-1] ->Fill(PhotEgy,PartThetaCM,yieldpol1);
		           fHistPiPlusAsym2Egphi[1][polstate] ->Fill(PhotEgy,PartPhi,yieldpol1);
		           fHistPiPlusAsym2Egtheta[1][polstate] ->Fill(PhotEgy,PartTheta,yieldpol1);
			   if (PhotEgy < MaxEgy)
			     {
		               fHistPiPlusAsym[1][polstate][tagChannel] ->Fill(PartPhi,PartTheta,yieldpol1);
			       fHistPiPlusAsymThetaface[1][polstate][face-1] ->Fill(PartTheta,yieldpol1);
			       fHistPiPlusAsym1Theta[1][polstate] ->Fill(PartTheta,yieldpol1);
			       fHistPiPlusAsym1Phi[1][polstate] ->Fill(PartPhi,yieldpol1);
		               fHistPiPlusAsym2Thetaphi[1][polstate] ->Fill(PartTheta,PartPhi,yieldpol1);
		               fHistPiPlusAsymHRThetaphi[1][polstate] ->Fill(PartTheta,PartPhi,yieldpol1);
		               fHistPiPlusAsymHRCosThetaphi[1][polstate] ->Fill(cos(RadPartTheta),PartPhi,yieldpol1);
			       fHistPiPlusAsym1ThetaCMface[1][polstate][face-1][tagChannel] ->Fill(PartThetaCM,yieldpol1);
			       fHistPiPlusAsym1Pion[1][face-1][tagChannel] ->Fill(partEgy,1);
			       fHistPiPlusAsym2PiTheta[1][face-1][tagChannel] ->Fill(partEgy,PartTheta,1);
			       fHistPiPlusAsym2PiPhi[1][tagChannel] ->Fill(partEgy,PartPhi,1);
			       fHistPiPlusAsym3PiThePhi[1][tagChannel] ->Fill(partEgy,PartTheta,PartPhi,1);
			     }

			   fHistPiPlusAsym1Eg[2][polstate] ->Fill(PhotEgy);
			   fHistPiPlusAsymEgface[2][polstate][face-1] ->Fill(PhotEgy);
		           fHistPiPlusAsym2EgthetaLBface[2][polstate][face-1] ->Fill(PhotEgy,PartTheta);
		           fHistPiPlusAsym2EgthetaCMface[2][polstate][face-1] ->Fill(PhotEgy,PartThetaCM);
		           fHistPiPlusAsym2Egphi[2][polstate] ->Fill(PhotEgy,PartPhi);
		           fHistPiPlusAsym2Egtheta[2][polstate] ->Fill(PhotEgy,PartTheta);
			   if (PhotEgy < MaxEgy)
			     {
		               fHistPiPlusAsym[2][polstate][tagChannel] ->Fill(PartPhi,PartTheta);
			       fHistPiPlusAsymThetaface[2][polstate][face-1] ->Fill(PartTheta);
			       fHistPiPlusAsym1Theta[2][polstate] ->Fill(PartTheta);
			       fHistPiPlusAsym1Phi[2][polstate] ->Fill(PartPhi);
		               fHistPiPlusAsym2Thetaphi[2][polstate] ->Fill(PartTheta,PartPhi);
		               fHistPiPlusAsymHRThetaphi[2][polstate] ->Fill(PartTheta,PartPhi);
		               fHistPiPlusAsymHRCosThetaphi[2][polstate] ->Fill(cos(RadPartTheta),PartPhi);
			       fHistPiPlusAsym1ThetaCMface[2][polstate][face-1][tagChannel] ->Fill(PartThetaCM);
			       fHistPiPlusAsym1Pion[2][face-1][tagChannel] ->Fill(partEgy);
			       fHistPiPlusAsym2PiTheta[2][face-1][tagChannel] ->Fill(partEgy,PartTheta);
			       fHistPiPlusAsym2PiPhi[2][tagChannel] ->Fill(partEgy,PartPhi);
			       fHistPiPlusAsym3PiThePhi[2][tagChannel] ->Fill(partEgy,PartTheta,PartPhi);
			     }
			}
		      else if (tagtdc >= HighTrueTag && tagtdc < HighAccTag)
			{
			   fHistPiPlusAsym1Eg[0][polstate] ->Fill(PhotEgy,-norm1acc);
			   fHistPiPlusAsymEgface[0][polstate][face-1] ->Fill(PhotEgy,-norm1acc);
		           fHistPiPlusAsym2EgthetaLBface[0][polstate][face-1] ->Fill(PhotEgy,PartTheta,-norm1acc);
		           fHistPiPlusAsym2EgthetaCMface[0][polstate][face-1] ->Fill(PhotEgy,PartThetaCM,-norm1acc);
		           fHistPiPlusAsym2Egphi[0][polstate] ->Fill(PhotEgy,PartPhi,-norm1acc);
		           fHistPiPlusAsym2Egtheta[0][polstate] ->Fill(PhotEgy,PartTheta,-norm1acc);
			   if (PhotEgy < MaxEgy)
			     {
		               fHistPiPlusAsym[0][polstate][tagChannel] ->Fill(PartPhi,PartTheta,-norm1acc);
			       fHistPiPlusAsymThetaface[0][polstate][face-1] ->Fill(PartTheta,-norm1acc);
			       fHistPiPlusAsym1Theta[0][polstate] ->Fill(PartTheta,-norm1acc);
			       fHistPiPlusAsym1Phi[0][polstate] ->Fill(PartPhi,-norm1acc);
		               fHistPiPlusAsym2Thetaphi[0][polstate] ->Fill(PartTheta,PartPhi,-norm1acc);
		               fHistPiPlusAsymHRThetaphi[0][polstate] ->Fill(PartTheta,PartPhi,-norm1acc);
		               fHistPiPlusAsymHRCosThetaphi[0][polstate] ->Fill(cos(RadPartTheta),PartPhi,-norm1acc);
			       fHistPiPlusAsym1ThetaCMface[0][polstate][face-1][tagChannel] ->Fill(PartThetaCM,-norm1acc);
			       fHistPiPlusAsym1Pion[0][face-1][tagChannel] ->Fill(partEgy,-accscal);
			       fHistPiPlusAsym2PiTheta[0][face-1][tagChannel] ->Fill(partEgy,PartTheta,-accscal);
			       fHistPiPlusAsym2PiPhi[0][tagChannel] ->Fill(partEgy,PartPhi,-accscal);
			       fHistPiPlusAsym3PiThePhi[0][tagChannel] ->Fill(partEgy,PartTheta,PartPhi,-accscal);
			     }

			   fHistPiPlusAsym1Eg[1][polstate] ->Fill(PhotEgy,-yieldpol1acc);
			   fHistPiPlusAsymEgface[1][polstate][face-1] ->Fill(PhotEgy,-yieldpol1acc);
		           fHistPiPlusAsym2EgthetaLBface[1][polstate][face-1] ->Fill(PhotEgy,PartTheta,-yieldpol1acc);
		           fHistPiPlusAsym2EgthetaCMface[1][polstate][face-1] ->Fill(PhotEgy,PartThetaCM,-yieldpol1acc);
		           fHistPiPlusAsym2Egphi[1][polstate] ->Fill(PhotEgy,PartPhi,-yieldpol1acc);
		           fHistPiPlusAsym2Egtheta[1][polstate] ->Fill(PhotEgy,PartTheta,-yieldpol1acc);
			   if (PhotEgy < MaxEgy)
			     {
		               fHistPiPlusAsym[1][polstate][tagChannel] ->Fill(PartPhi,PartTheta,-yieldpol1acc);
			       fHistPiPlusAsymThetaface[1][polstate][face-1] ->Fill(PartTheta,-yieldpol1acc);
			       fHistPiPlusAsym1Theta[1][polstate] ->Fill(PartTheta,-yieldpol1acc);
			       fHistPiPlusAsym1Phi[1][polstate] ->Fill(PartPhi,-yieldpol1acc);
		               fHistPiPlusAsym2Thetaphi[1][polstate] ->Fill(PartTheta,PartPhi,-yieldpol1acc);
		               fHistPiPlusAsymHRThetaphi[1][polstate] ->Fill(PartTheta,PartPhi,-yieldpol1acc);
		               fHistPiPlusAsymHRCosThetaphi[1][polstate] ->Fill(cos(RadPartTheta),PartPhi,-yieldpol1acc);
			       fHistPiPlusAsym1ThetaCMface[1][polstate][face-1][tagChannel] ->Fill(PartThetaCM,-yieldpol1acc);
			       fHistPiPlusAsym1Pion[1][face-1][tagChannel] ->Fill(partEgy,-accscal);
			       fHistPiPlusAsym2PiTheta[1][face-1][tagChannel] ->Fill(partEgy,PartTheta,-accscal);
			       fHistPiPlusAsym2PiPhi[1][tagChannel] ->Fill(partEgy,PartPhi,-accscal);
			       fHistPiPlusAsym3PiThePhi[1][tagChannel] ->Fill(partEgy,PartTheta,PartPhi,-accscal);
			     }

			   fHistPiPlusAsym1Eg[3][polstate] ->Fill(PhotEgy);
			   fHistPiPlusAsymEgface[3][polstate][face-1] ->Fill(PhotEgy);
		           fHistPiPlusAsym2EgthetaLBface[3][polstate][face-1] ->Fill(PhotEgy,PartTheta);
		           fHistPiPlusAsym2EgthetaCMface[3][polstate][face-1] ->Fill(PhotEgy,PartThetaCM);
		           fHistPiPlusAsym2Egphi[3][polstate] ->Fill(PhotEgy,PartPhi);
		           fHistPiPlusAsym2Egtheta[3][polstate] ->Fill(PhotEgy,PartTheta);
			   if (PhotEgy < MaxEgy)
			     {
		               fHistPiPlusAsym[3][polstate][tagChannel] ->Fill(PartPhi,PartTheta);
			       fHistPiPlusAsymThetaface[3][polstate][face-1] ->Fill(PartTheta);
			       fHistPiPlusAsym1Theta[3][polstate] ->Fill(PartTheta);
			       fHistPiPlusAsym1Phi[3][polstate] ->Fill(PartPhi);
		               fHistPiPlusAsym2Thetaphi[3][polstate] ->Fill(PartTheta,PartPhi);
		               fHistPiPlusAsymHRThetaphi[3][polstate] ->Fill(PartTheta,PartPhi);
		               fHistPiPlusAsymHRCosThetaphi[3][polstate] ->Fill(cos(RadPartTheta),PartPhi);
			       fHistPiPlusAsym1ThetaCMface[3][polstate][face-1][tagChannel] ->Fill(PartThetaCM);
			       fHistPiPlusAsym1Pion[3][face-1][tagChannel] ->Fill(partEgy);
			       fHistPiPlusAsym2PiTheta[3][face-1][tagChannel] ->Fill(partEgy,PartTheta);
			       fHistPiPlusAsym2PiPhi[3][tagChannel] ->Fill(partEgy,PartPhi);
			       fHistPiPlusAsym3PiThePhi[3][tagChannel] ->Fill(partEgy,PartTheta,PartPhi);
			     }
			}
		      else
			{
			   fHistPiPlusAsym1Eg[4][polstate] ->Fill(PhotEgy);
			   fHistPiPlusAsymEgface[4][polstate][face-1] ->Fill(PhotEgy);
		           fHistPiPlusAsym2EgthetaLBface[4][polstate][face-1] ->Fill(PhotEgy,PartTheta);
		           fHistPiPlusAsym2EgthetaCMface[4][polstate][face-1] ->Fill(PhotEgy,PartThetaCM);
		           fHistPiPlusAsym2Egphi[4][polstate] ->Fill(PhotEgy,PartPhi);
		           fHistPiPlusAsym2Egtheta[4][polstate] ->Fill(PhotEgy,PartTheta);
			   if (PhotEgy < MaxEgy)
			     {
		               fHistPiPlusAsym[4][polstate][tagChannel] ->Fill(PartPhi,PartTheta);
			       fHistPiPlusAsymThetaface[4][polstate][face-1] ->Fill(PartTheta);
			       fHistPiPlusAsym1Theta[4][polstate] ->Fill(PartTheta);
			       fHistPiPlusAsym1Phi[4][polstate] ->Fill(PartPhi);
		               fHistPiPlusAsym2Thetaphi[4][polstate] ->Fill(PartTheta,PartPhi);
		               fHistPiPlusAsymHRThetaphi[4][polstate] ->Fill(PartTheta,PartPhi);
		               fHistPiPlusAsymHRCosThetaphi[4][polstate] ->Fill(cos(RadPartTheta),PartPhi);
			       fHistPiPlusAsym1ThetaCMface[4][polstate][face-1][tagChannel] ->Fill(PartTheta);
			       fHistPiPlusAsym1Pion[4][face-1][tagChannel] ->Fill(partEgy);
			       fHistPiPlusAsym2PiTheta[4][face-1][tagChannel] ->Fill(partEgy,PartTheta);
			       fHistPiPlusAsym2PiPhi[4][tagChannel] ->Fill(partEgy,PartPhi);
			       fHistPiPlusAsym3PiThePhi[4][tagChannel] ->Fill(partEgy,PartTheta,PartPhi);
			     }
			}
		}
	    }
	}
    }
  return flag;
}



ROOT page - Home page - Class index - Class Hierarchy - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.