1390c238672c6fdda6b98589c66cc23ea4c77377
[u/mrichter/AliRoot.git] / MUON / MUONTestAbso.C
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(pz0*pz0+px0*px0);
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(pYZ*pYZ+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(x*x+y*y);
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 }