[u/mrichter/AliRoot.git] / MUON / MUONTestAbso.C

void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0) 
{
  // Useful to test the absorber correction in the reconstruction procedure
  // (mean energy loss + Branson correction)
  // The AliMUONTrackParam class from the reconstruction is directly checked
  // in this macro using the GEANT particle momentum downstream of the absorber.
  // Histograms are saved in file MUONTestAbso.root.
  // Use DrawTestAbso.C to display control histograms.

    
  // Dynamically link some shared libs                    
  if (gClassTable->GetID("AliRun") < 0) {
    gROOT->LoadMacro("loadlibs.C");
    loadlibs();
  }

  // Connect the Root Galice file containing Geometry, Kine and Hits

  TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");

    
  if (!file) {
    printf("\n Creating galice.root \n");
    file = new TFile("galice.root");
  } else {
    printf("\n galice.root found in file list");
  }
    
// Get AliRun object from file or create it if not on file
  if (!gAlice) {
    gAlice = (AliRun*)(file->Get("gAlice"));
    if (gAlice) printf("AliRun object found on file\n");
    if (!gAlice) {
      printf("\n Create new gAlice object");
      gAlice = new AliRun("gAlice","Alice test program");
    }
  }

    
  //  Create some histograms

  TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 120, 8., 11.);
  TH1F *hDeltaP1 = new TH1F("hDeltaP1", " delta P for muons r < rLimit", 100, -10., 10.);
  TH1F *hDeltaAngle1 = new TH1F("hDeltaAngle1", " delta angle for muons r < rLimit", 100, 0., 0.5);
  //
  TH1F *hDeltaP2 = new TH1F("hDeltaP2", " delta P for muons r > rLimit", 100, -10., 10.);
  TH1F *hDeltaAngle2 = new TH1F("hDeltaAngle2", " delta angle for muons r > rLimit", 100, 0., 0.5);
  
  // Define constant parameters
  Float_t massMuon = 0.105658389;  // muon mass
  Float_t massUpsilon = 9.46037;   // Upsilon mass
  Double_t PI = 3.141592654;
  Double_t zEndAbso = 503;         // z position of the end of the absorber 
  Double_t rLimit = zEndAbso * TMath::Tan(3*PI/180); // 3 deg. limit (different absorber composition)
  Double_t printLevel = 0;

  if (printLevel >= 1)
  cout <<" **** z end Absorber : "<< zEndAbso <<" rLimit : "<< rLimit << endl;

  Double_t pX[2],pY[2],pZ[2],pTot[2],theta[2],radius[2]; // extrapolated parameters of particles from chamber 1 to vertex 
   
  for (Int_t nev=0; nev<= evNumber2; nev++) {  // loop over events
    if (printLevel >= 0) {
      cout << "  "<<endl;
      cout << " **** Event # " << nev <<endl;
    }
    
    Int_t nparticles = gAlice->GetEvent(nev);
    if (printLevel >= 1) {
      cout << " Total number of nparticles  " << nparticles <<endl;
    }
    if (nev < evNumber1) continue;
    if (nparticles <= 0) return;
     
    Double_t PxGen[2],PyGen[2],PzGen[2],PTotGen[2],ThetaGen[2]; // parameters of generated  particles at the vertex
     
    for (int iPart = 0; iPart < nparticles  ; iPart++) {  // loop over particles
      
      TParticle *particle;
      particle = gAlice->Particle(iPart);    

      if ( particle->GetPdgCode() == kMuonPlus ) {
	PxGen[0] = particle->Px();
	PyGen[0] = particle->Py();
	PzGen[0] = particle->Pz();
	PTotGen[0] = TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]+PzGen[0]*PzGen[0]);
	ThetaGen[0] = TMath::ATan2(TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]),PzGen[0])*180/PI;
	if (printLevel >= 1) {
	  cout << " Particle id : "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[0]<<" "<<PyGen[0]<<" "<<PzGen[0]<<" "<<ThetaGen[0]<<endl;
	}
      }
      if ( particle->GetPdgCode() == kMuonMinus ) {
	PxGen[1] = particle->Px();
	PyGen[1] = particle->Py();
	PzGen[1] = particle->Pz();
	PTotGen[1] = TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]+PzGen[1]*PzGen[1]);
	ThetaGen[1] = TMath::ATan2(TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]),PzGen[1])*180/PI;
	if (printLevel >= 1)
	  cout << " Particle #: "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[1]<<" "<<PyGen[1]<<" "<<PzGen[1]<<" "<<ThetaGen[1]<<endl;
      }
    }


    AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
     
     
    TTree *TH = gAlice->TreeH();
    Int_t ntracks =(Int_t) TH->GetEntries();
     
    for (Int_t i = 0; i < 2; i++){  
      pX[i] = pY[i] = pZ[i] = pTot[i] = theta[i] = radius[i] = 0;
    }
    
    for (Int_t track=0; track<ntracks;track++) { // loop over tracks
      gAlice->ResetHits();
      Int_t nbytes += TH->GetEvent(track);
      if (MUON)  {
	for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); 
	    mHit;
	    mHit=(AliMUONHit*)MUON->NextHit())  // loop over GEANT muon hits
	  {
	    Int_t   nch   = mHit->Chamber();  // chamber number
	     
	    if (nch != 1) continue;  
	    
	    Int_t ipart = mHit->Particle();
	    Double_t x     = mHit->X();        // x-pos of hit in chamber #1
	    Double_t y     = mHit->Y();        // y-pos of hit in chamber #1
	    Double_t z     = mHit->Z();        // z-pos of hit in chamber #1
	     
	     
	    if (ipart == kMuonPlus || ipart == kMuonMinus) {
	      Double_t px0=mHit->Px();
	      Double_t py0=mHit->Py();		   
	      Double_t pz0=mHit->Pz();
	      Double_t nonBendingSlope=px0/pz0;
	      Double_t bendingSlope=py0/pz0;

	      // create an AliMUONTrackParam object in chamber #1
	      AliMUONTrackParam *trackParam = new AliMUONTrackParam();
	      Double_t bendingMometum = TMath::Sqrt(pz0*pz0+px0*px0);
	      Double_t inverseBendingMomentum = 1/bendingMometum;
	       
	      trackParam->SetBendingCoor(y);
	      trackParam->SetNonBendingCoor(x);
	      trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
	      trackParam->SetBendingSlope(bendingSlope);
	      trackParam->SetNonBendingSlope(nonBendingSlope);
	      trackParam->SetZ(z);	
	       
	      if (printLevel >= 1) {
		cout <<" **** before extrap " <<endl;
		cout << " px, py, pz = " <<px0<<" "<<py0<<" "<<pz0<<endl;
		trackParam->Dump();
	      }
	      trackParam->ExtrapToVertex(); // extrapolate track parameters through the absorber
	      if (printLevel >= 1) {
		cout <<" **** after extrap " <<endl;
		trackParam->Dump();
	      }

	      Int_t iMuon;
	      if (ipart == kMuonPlus) iMuon = 0;
	      else iMuon = 1;
	       
	      // calculate track parameters extrapolated to vertex.
	      bendingSlope = trackParam->GetBendingSlope();
	      nonBendingSlope = trackParam->GetNonBendingSlope();
	      Double_t pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
	      pZ[iMuon] = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
	      pX[iMuon] = pZ[iMuon] * nonBendingSlope;
	      pY[iMuon] = pZ[iMuon] * bendingSlope;
	      pTot[iMuon] = TMath::Sqrt(pYZ*pYZ+pX[iMuon]*pX[iMuon]);
	      theta[iMuon] = TMath::ATan2(TMath::Sqrt(pX[iMuon]*pX[iMuon]+pY[iMuon]*pY[iMuon]),pZ[iMuon])*180/PI;
	      radius[iMuon] = TMath::Sqrt(x*x+y*y);
	      
	      if (printLevel >= 1)
		cout << " px, py, pz, theta, radius = " <<pX[iMuon]<<" "<<pY[iMuon]<<" "<<pZ[iMuon]<<" "<<theta[iMuon]<<" "<<radius[iMuon]<<endl;
	      delete trackParam;
	      
	    } // if MuonPlus or MuonMinus
	  } // loop over MUON hits 
      } // if MUON 
    } // loop over tracks	   

    if (pX[0] != 0 && pX[1] != 0) { // if track 1 & 2 exist 
      Float_t massResonance = MuPlusMuMinusMass(pX[0],pY[0],pZ[0],pX[1],pY[1],pZ[1],massMuon);
      hInvMassRes->Fill(massResonance);
      for (Int_t i = 0; i<2; i++){
	Double_t cosAngle = pX[i]*PxGen[i]+pY[i]*PyGen[i]+pZ[i]*PzGen[i];
	cosAngle = cosAngle/(pTot[i]*PTotGen[i]);
	Double_t deltaAngle = TMath::ACos(cosAngle)*180/PI;
	if (radius[i] < rLimit) {
	  hDeltaP1->Fill(pTot[i]-PTotGen[i]);
	  hDeltaAngle1->Fill(deltaAngle);
	}else{
	  hDeltaP2->Fill(pTot[i]-PTotGen[i]);
	  hDeltaAngle2->Fill(deltaAngle);
	}
      }
    }
  } // loop over event

  // save histograms in MUONTestAbso.root
  TFile *histoFile = new TFile("MUONTestAbso.root", "RECREATE");
  hInvMassRes->Write();
  hDeltaP1->Write();
  hDeltaAngle1->Write();
  hDeltaP2->Write();
  hDeltaAngle2->Write();
  histoFile->Close();
  
}


Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t MuonMass)
{
  // return invariant mass for particle 1 & 2 whose masses are equal to MuonMass

  Float_t e1 = TMath::Sqrt(MuonMass * MuonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
  Float_t e2 = TMath::Sqrt(MuonMass * MuonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
  return (TMath::Sqrt(2.0 * (MuonMass * MuonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
}
Commit	Line	Data
b6340d07	1	void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
	2	{
	3	// Useful to test the absorber correction in the reconstruction procedure
	4	// (mean energy loss + Branson correction)
	5	// The AliMUONTrackParam class from the reconstruction is directly checked
	6	// in this macro using the GEANT particle momentum downstream of the absorber.
	7	// Histograms are saved in file MUONTestAbso.root.
	8	// Use DrawTestAbso.C to display control histograms.
	9
	10
	11	// Dynamically link some shared libs
	12	if (gClassTable->GetID("AliRun") < 0) {
	13	gROOT->LoadMacro("loadlibs.C");
	14	loadlibs();
	15	}
	16
	17	// Connect the Root Galice file containing Geometry, Kine and Hits
	18
	19	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject("galice.root");
	20
	21
	22	if (!file) {
	23	printf("\n Creating galice.root \n");
	24	file = new TFile("galice.root");
	25	} else {
	26	printf("\n galice.root found in file list");
	27	}
	28
	29	// Get AliRun object from file or create it if not on file
	30	if (!gAlice) {
	31	gAlice = (AliRun*)(file->Get("gAlice"));
	32	if (gAlice) printf("AliRun object found on file\n");
	33	if (!gAlice) {
	34	printf("\n Create new gAlice object");
	35	gAlice = new AliRun("gAlice","Alice test program");
	36	}
	37	}
	38
	39
	40	// Create some histograms
	41
	42	TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 120, 8., 11.);
	43	TH1F *hDeltaP1 = new TH1F("hDeltaP1", " delta P for muons r < rLimit", 100, -10., 10.);
	44	TH1F *hDeltaAngle1 = new TH1F("hDeltaAngle1", " delta angle for muons r < rLimit", 100, 0., 0.5);
	45	//
	46	TH1F *hDeltaP2 = new TH1F("hDeltaP2", " delta P for muons r > rLimit", 100, -10., 10.);
	47	TH1F *hDeltaAngle2 = new TH1F("hDeltaAngle2", " delta angle for muons r > rLimit", 100, 0., 0.5);
	48
	49	// Define constant parameters
	50	Float_t massMuon = 0.105658389; // muon mass
	51	Float_t massUpsilon = 9.46037; // Upsilon mass
	52	Double_t PI = 3.141592654;
	53	Double_t zEndAbso = 503; // z position of the end of the absorber
	54	Double_t rLimit = zEndAbso * TMath::Tan(3*PI/180); // 3 deg. limit (different absorber composition)
	55	Double_t printLevel = 0;
	56
	57	if (printLevel >= 1)
	58	cout <<" **** z end Absorber : "<< zEndAbso <<" rLimit : "<< rLimit << endl;
	59
	60	Double_t pX[2],pY[2],pZ[2],pTot[2],theta[2],radius[2]; // extrapolated parameters of particles from chamber 1 to vertex
	61
	62	for (Int_t nev=0; nev<= evNumber2; nev++) { // loop over events
	63	if (printLevel >= 0) {
	64	cout << " "<<endl;
65	cout << " **** Event # " << nev <<endl;
66	}
67
68	Int_t nparticles = gAlice->GetEvent(nev);
69	if (printLevel >= 1) {
70	cout << " Total number of nparticles " << nparticles <<endl;
71	}
72	if (nev < evNumber1) continue;
73	if (nparticles <= 0) return;
74
75	Double_t PxGen[2],PyGen[2],PzGen[2],PTotGen[2],ThetaGen[2]; // parameters of generated particles at the vertex
76
77	for (int iPart = 0; iPart < nparticles ; iPart++) { // loop over particles
78
79	TParticle *particle;
80	particle = gAlice->Particle(iPart);
81
82	if ( particle->GetPdgCode() == kMuonPlus ) {
83	PxGen[0] = particle->Px();
84	PyGen[0] = particle->Py();
85	PzGen[0] = particle->Pz();
86	PTotGen[0] = TMath::Sqrt(PxGen[0]PxGen[0]+PyGen[0]PyGen[0]+PzGen[0]*PzGen[0]);
87	ThetaGen[0] = TMath::ATan2(TMath::Sqrt(PxGen[0]PxGen[0]+PyGen[0]PyGen[0]),PzGen[0])*180/PI;
88	if (printLevel >= 1) {
89	cout << " Particle id : "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[0]<<" "<<PyGen[0]<<" "<<PzGen[0]<<" "<<ThetaGen[0]<<endl;
90	}
91	}
92	if ( particle->GetPdgCode() == kMuonMinus ) {
93	PxGen[1] = particle->Px();
94	PyGen[1] = particle->Py();
95	PzGen[1] = particle->Pz();
96	PTotGen[1] = TMath::Sqrt(PxGen[1]PxGen[1]+PyGen[1]PyGen[1]+PzGen[1]*PzGen[1]);
97	ThetaGen[1] = TMath::ATan2(TMath::Sqrt(PxGen[1]PxGen[1]+PyGen[1]PyGen[1]),PzGen[1])*180/PI;
98	if (printLevel >= 1)
99	cout << " Particle #: "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[1]<<" "<<PyGen[1]<<" "<<PzGen[1]<<" "<<ThetaGen[1]<<endl;
100	}
101	}
102
103
104	AliMUON MUON = (AliMUON) gAlice->GetModule("MUON");
105
106
107	TTree *TH = gAlice->TreeH();
108	Int_t ntracks =(Int_t) TH->GetEntries();
109
110	for (Int_t i = 0; i < 2; i++){
111	pX[i] = pY[i] = pZ[i] = pTot[i] = theta[i] = radius[i] = 0;
112	}
113
114	for (Int_t track=0; track<ntracks;track++) { // loop over tracks
115	gAlice->ResetHits();
116	Int_t nbytes += TH->GetEvent(track);
117	if (MUON) {
118	for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
119	mHit;
120	mHit=(AliMUONHit*)MUON->NextHit()) // loop over GEANT muon hits
121	{
122	Int_t nch = mHit->Chamber(); // chamber number
123
124	if (nch != 1) continue;
125
126	Int_t ipart = mHit->Particle();
127	Double_t x = mHit->X(); // x-pos of hit in chamber #1
128	Double_t y = mHit->Y(); // y-pos of hit in chamber #1
129	Double_t z = mHit->Z(); // z-pos of hit in chamber #1
130
131
132	if (ipart == kMuonPlus \|\| ipart == kMuonMinus) {
133	Double_t px0=mHit->Px();
134	Double_t py0=mHit->Py();
135	Double_t pz0=mHit->Pz();
136	Double_t nonBendingSlope=px0/pz0;
137	Double_t bendingSlope=py0/pz0;
138
139	// create an AliMUONTrackParam object in chamber #1
140	AliMUONTrackParam *trackParam = new AliMUONTrackParam();
141	Double_t bendingMometum = TMath::Sqrt(pz0pz0+px0px0);
142	Double_t inverseBendingMomentum = 1/bendingMometum;
143
144	trackParam->SetBendingCoor(y);
145	trackParam->SetNonBendingCoor(x);
146	trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
147	trackParam->SetBendingSlope(bendingSlope);
148	trackParam->SetNonBendingSlope(nonBendingSlope);
149	trackParam->SetZ(z);
150
151	if (printLevel >= 1) {
152	cout <<" **** before extrap " <<endl;
153	cout << " px, py, pz = " <<px0<<" "<<py0<<" "<<pz0<<endl;
154	trackParam->Dump();
155	}
156	trackParam->ExtrapToVertex(); // extrapolate track parameters through the absorber
157	if (printLevel >= 1) {
158	cout <<" **** after extrap " <<endl;
159	trackParam->Dump();
160	}
161
162	Int_t iMuon;
163	if (ipart == kMuonPlus) iMuon = 0;
164	else iMuon = 1;
165
166	// calculate track parameters extrapolated to vertex.
167	bendingSlope = trackParam->GetBendingSlope();
168	nonBendingSlope = trackParam->GetNonBendingSlope();
169	Double_t pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
170	pZ[iMuon] = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
171	pX[iMuon] = pZ[iMuon] * nonBendingSlope;
172	pY[iMuon] = pZ[iMuon] * bendingSlope;
173	pTot[iMuon] = TMath::Sqrt(pYZpYZ+pX[iMuon]pX[iMuon]);
174	theta[iMuon] = TMath::ATan2(TMath::Sqrt(pX[iMuon]pX[iMuon]+pY[iMuon]pY[iMuon]),pZ[iMuon])*180/PI;
175	radius[iMuon] = TMath::Sqrt(xx+yy);
176
177	if (printLevel >= 1)
178	cout << " px, py, pz, theta, radius = " <<pX[iMuon]<<" "<<pY[iMuon]<<" "<<pZ[iMuon]<<" "<<theta[iMuon]<<" "<<radius[iMuon]<<endl;
179	delete trackParam;
180
181	} // if MuonPlus or MuonMinus
182	} // loop over MUON hits
183	} // if MUON
184	} // loop over tracks
185
186	if (pX[0] != 0 && pX[1] != 0) { // if track 1 & 2 exist
187	Float_t massResonance = MuPlusMuMinusMass(pX[0],pY[0],pZ[0],pX[1],pY[1],pZ[1],massMuon);
188	hInvMassRes->Fill(massResonance);
189	for (Int_t i = 0; i<2; i++){
190	Double_t cosAngle = pX[i]PxGen[i]+pY[i]PyGen[i]+pZ[i]*PzGen[i];
191	cosAngle = cosAngle/(pTot[i]*PTotGen[i]);
192	Double_t deltaAngle = TMath::ACos(cosAngle)*180/PI;
193	if (radius[i] < rLimit) {
194	hDeltaP1->Fill(pTot[i]-PTotGen[i]);
195	hDeltaAngle1->Fill(deltaAngle);
196	}else{
197	hDeltaP2->Fill(pTot[i]-PTotGen[i]);
198	hDeltaAngle2->Fill(deltaAngle);
199	}
200	}
201	}
202	} // loop over event
203
204	// save histograms in MUONTestAbso.root
205	TFile *histoFile = new TFile("MUONTestAbso.root", "RECREATE");
206	hInvMassRes->Write();
207	hDeltaP1->Write();
208	hDeltaAngle1->Write();
209	hDeltaP2->Write();
210	hDeltaAngle2->Write();
211	histoFile->Close();
212
213	}
214
215
216	Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t MuonMass)
217	{
218	// return invariant mass for particle 1 & 2 whose masses are equal to MuonMass
219
220	Float_t e1 = TMath::Sqrt(MuonMass * MuonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
221	Float_t e2 = TMath::Sqrt(MuonMass * MuonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
222	return (TMath::Sqrt(2.0 * (MuonMass * MuonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
223	}