Write the parameter class into the cluster file
[u/mrichter/AliRoot.git] / MUON / MUONTestAbso.C
CommitLineData
b6340d07 1void 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
216Float_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}