]>
Commit | Line | Data |
---|---|---|
b45fd22b | 1 | void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0,Int_t testSingle=0) |
b6340d07 | 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 | ||
b45fd22b | 42 | TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2)", 1200, 0., 12.); |
43 | // TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2)", 100, 9., 10.5); | |
44 | TH1F *hDeltaP1 = new TH1F("hDeltaP1", " delta P for muons theta < 3 deg", 100, -10., 10.); | |
45 | TH1F *hDeltaAngle1 = new TH1F("hDeltaAngle1", " delta angle for muons theta < 3 deg", 100, 0., 0.5); | |
b6340d07 | 46 | // |
b45fd22b | 47 | TH1F *hDeltaP2 = new TH1F("hDeltaP2", " delta P for muons theta > 3 deg", 100, -10., 10.); |
48 | TH1F *hDeltaAngle2 = new TH1F("hDeltaAngle2", " delta angle for muons theta > 3 deg", 100, 0., 0.5); | |
49 | ||
50 | // | |
51 | TH1F *hRap = new TH1F("hRap"," Muon Rapidity gen.",20,2.5,4); | |
52 | // | |
53 | ||
54 | Int_t nBinProf=19; | |
55 | Float_t xProfMin=5; | |
56 | Float_t xProfMax=290; | |
57 | char titleHist[50]; | |
58 | char numberHist[10]; | |
59 | TH1F *DeltaP1[100],*DeltaP2[100]; | |
60 | TH1F *DeltaAngleX1[100],*DeltaAngleX2[100]; | |
61 | TH1F *DeltaAngleY1[100],*DeltaAngleY2[100]; | |
b6340d07 | 62 | |
b45fd22b | 63 | TH1F *hP1 = new TH1F("hP1"," P theta < 3 deg ",nBinProf,xProfMin,xProfMax); |
64 | TProfile *hProfDeltaPvsP1 = new TProfile("hProfDeltaPvsP1"," delta P vs P theta < 3 deg ",nBinProf,xProfMin,xProfMax,-50,50,"s"); | |
65 | TH2F *h2DeltaPvsP1 = new TH2F("h2DeltaPvsP1"," delta P vs P theta < 3 deg",nBinProf,xProfMin,xProfMax,50,-50,50); | |
66 | TH1F *hSigmaPvsP1 = new TH1F("hSigmaPvsP1"," deltaP/P vs P theta < 3 deg",nBinProf,xProfMin,xProfMax); | |
67 | for (Int_t iHist=0; iHist < nBinProf; iHist++){ | |
68 | sprintf(titleHist,"%s","hDelta1P"); | |
69 | sprintf(numberHist,"%d",iHist+1); | |
70 | strcat(titleHist,numberHist); | |
71 | DeltaP1[iHist] = new TH1F(titleHist," deltaP theta < 3 deg ",400,-50,50); | |
72 | ||
73 | sprintf(titleHist,"%s","hDelta1AngleX"); | |
74 | sprintf(numberHist,"%d",iHist+1); | |
75 | strcat(titleHist,numberHist); | |
76 | DeltaAngleX1[iHist] = new TH1F(titleHist," delta Angle X theta < 3 deg ",400,-0.05,0.05); | |
77 | ||
78 | sprintf(titleHist,"%s","hDelta1AngleY"); | |
79 | sprintf(numberHist,"%d",iHist+1); | |
80 | strcat(titleHist,numberHist); | |
81 | DeltaAngleY1[iHist] = new TH1F(titleHist," delta Angle Y theta < 3 deg ",400,-0.05,0.05); | |
82 | } | |
83 | TH1F *hSigGausPvsP1 = new TH1F("hSigGausPvsP1"," deltaP/P vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax); | |
84 | TH1F *hSigGausAngleXvsP1 = new TH1F("hSigGausAngleXvsP1"," delta theta X vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax); | |
85 | TH1F *hSigGausAngleYvsP1 = new TH1F("hSigGausAngleYvsP1"," delta theta Y vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax); | |
86 | ||
87 | TH1F *hP2 = new TH1F("hP2"," P theta > 3 deg ",nBinProf,xProfMin,xProfMax); | |
88 | TProfile *hProfDeltaPvsP2 = new TProfile("hProfDeltaPvsP2"," delta P vs P theta > 3 deg",nBinProf,xProfMin,xProfMax,-50,50,"s"); | |
89 | TH2F *h2DeltaPvsP2 = new TH2F("h2DeltaPvsP2"," delta P vs P theta > 3 deg",nBinProf,xProfMin,xProfMax,50,-50,50); | |
90 | TH1F *hSigmaPvsP2 = new TH1F("hSigmaPvsP2"," deltaP/P vs P theta > 3 deg",nBinProf,xProfMin,xProfMax); | |
91 | for (Int_t iHist=0; iHist < nBinProf; iHist++){ | |
92 | sprintf(titleHist,"%s","hDelta2P"); | |
93 | sprintf(numberHist,"%d",iHist+1); | |
94 | strcat(titleHist,numberHist); | |
95 | DeltaP2[iHist] = new TH1F(titleHist," deltaP theta > 3 deg ",400,-50,50); | |
96 | ||
97 | sprintf(titleHist,"%s","hDelta2AngleX"); | |
98 | sprintf(numberHist,"%d",iHist+1); | |
99 | strcat(titleHist,numberHist); | |
100 | DeltaAngleX2[iHist] = new TH1F(titleHist," delta Angle X theta > 3 deg ",400,-0.05,0.05); | |
101 | ||
102 | sprintf(titleHist,"%s","hDelta2AngleY"); | |
103 | sprintf(numberHist,"%d",iHist+1); | |
104 | strcat(titleHist,numberHist); | |
105 | DeltaAngleY2[iHist] = new TH1F(titleHist," delta Angle Y theta > 3 deg ",400,-0.05,0.05); | |
106 | } | |
107 | TH1F *hSigGausPvsP2 = new TH1F("hSigGausPvsP2"," deltaP/P vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax); | |
108 | ||
109 | TH1F *hSigGausAngleXvsP2 = new TH1F("hSigGausAngleXvsP2"," delta theta X vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax); | |
110 | TH1F *hSigGausAngleYvsP2 = new TH1F("hSigGausAngleYvsP2"," delta theta Y vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax); | |
111 | ||
112 | TProfile *hDeltaPxvsPhi = new TProfile("hDeltaPxvsPhi"," delta Px vs Phi",45,-180,180,-1,1,"s"); | |
113 | TProfile *hDeltaPyvsPhi = new TProfile("hDeltaPyvsPhi"," delta Py vs Phi",45,-180,180,-1,1,"s"); | |
114 | TProfile *hDeltaPhivsPz = new TProfile("hDeltaPhivsPz"," delta phi vs pZ",25,0,100,-4,4); | |
115 | ||
b6340d07 | 116 | // Define constant parameters |
117 | Float_t massMuon = 0.105658389; // muon mass | |
118 | Float_t massUpsilon = 9.46037; // Upsilon mass | |
b6340d07 | 119 | Double_t zEndAbso = 503; // z position of the end of the absorber |
b45fd22b | 120 | Double_t rLimit = zEndAbso * TMath::Tan(3*TMath::Pi()/180); // 3 deg. limit (different absorber composition) |
b6340d07 | 121 | Double_t printLevel = 0; |
122 | ||
123 | if (printLevel >= 1) | |
124 | cout <<" **** z end Absorber : "<< zEndAbso <<" rLimit : "<< rLimit << endl; | |
125 | ||
126 | Double_t pX[2],pY[2],pZ[2],pTot[2],theta[2],radius[2]; // extrapolated parameters of particles from chamber 1 to vertex | |
127 | ||
b45fd22b | 128 | for (Int_t nev=evNumber1; nev<= evNumber2; nev++) { // loop over events |
129 | if (printLevel >= 0 && nev%100 == 0) { | |
b6340d07 | 130 | cout << " **** Event # " << nev <<endl; |
131 | } | |
132 | ||
133 | Int_t nparticles = gAlice->GetEvent(nev); | |
134 | if (printLevel >= 1) { | |
135 | cout << " Total number of nparticles " << nparticles <<endl; | |
136 | } | |
137 | if (nev < evNumber1) continue; | |
138 | if (nparticles <= 0) return; | |
139 | ||
b45fd22b | 140 | Double_t PxGen[2],PyGen[2],PzGen[2],PTotGen[2],ThetaGen[2],RapGen[2]; // parameters of generated particles at the vertex |
141 | ThetaGen[0] = ThetaGen[1] = 0; | |
142 | RapGen[0] = RapGen[1] = 0; | |
143 | ||
b6340d07 | 144 | for (int iPart = 0; iPart < nparticles ; iPart++) { // loop over particles |
145 | ||
146 | TParticle *particle; | |
147 | particle = gAlice->Particle(iPart); | |
148 | ||
149 | if ( particle->GetPdgCode() == kMuonPlus ) { | |
150 | PxGen[0] = particle->Px(); | |
151 | PyGen[0] = particle->Py(); | |
152 | PzGen[0] = particle->Pz(); | |
153 | PTotGen[0] = TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]+PzGen[0]*PzGen[0]); | |
b45fd22b | 154 | ThetaGen[0] = TMath::ATan2(TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]),PzGen[0]); |
155 | RapGen[0] = rapParticle(PxGen[0],PyGen[0],PzGen[0],massMuon); | |
b6340d07 | 156 | if (printLevel >= 1) { |
157 | cout << " Particle id : "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[0]<<" "<<PyGen[0]<<" "<<PzGen[0]<<" "<<ThetaGen[0]<<endl; | |
158 | } | |
159 | } | |
160 | if ( particle->GetPdgCode() == kMuonMinus ) { | |
161 | PxGen[1] = particle->Px(); | |
162 | PyGen[1] = particle->Py(); | |
163 | PzGen[1] = particle->Pz(); | |
164 | PTotGen[1] = TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]+PzGen[1]*PzGen[1]); | |
b45fd22b | 165 | ThetaGen[1] = TMath::ATan2(TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]),PzGen[1]); |
166 | RapGen[1] = rapParticle(PxGen[1],PyGen[1],PzGen[1],massMuon); | |
b6340d07 | 167 | if (printLevel >= 1) |
168 | cout << " Particle #: "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[1]<<" "<<PyGen[1]<<" "<<PzGen[1]<<" "<<ThetaGen[1]<<endl; | |
169 | } | |
b45fd22b | 170 | } // end loop over particles |
b6340d07 | 171 | |
172 | ||
173 | AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); | |
174 | ||
175 | ||
176 | TTree *TH = gAlice->TreeH(); | |
177 | Int_t ntracks =(Int_t) TH->GetEntries(); | |
b45fd22b | 178 | |
179 | Bool_t hitChamberMUPlus[10],hitChamberMUMinus[10]; | |
180 | Bool_t hitStationMUPlus[5],hitStationMUMinus[5]; | |
181 | for (Int_t iCh=0; iCh < 10; iCh++){ | |
182 | hitChamberMUPlus[iCh] = kFALSE; | |
183 | hitChamberMUMinus[iCh] = kFALSE; | |
184 | } | |
185 | for (Int_t iSt=0; iSt < 5; iSt++){ | |
186 | hitStationMUPlus[iSt] = kFALSE; | |
187 | hitStationMUMinus[iSt] = kFALSE; | |
188 | } | |
189 | ||
190 | Int_t nHitChamber=0; | |
191 | ||
b6340d07 | 192 | for (Int_t i = 0; i < 2; i++){ |
193 | pX[i] = pY[i] = pZ[i] = pTot[i] = theta[i] = radius[i] = 0; | |
194 | } | |
195 | ||
196 | for (Int_t track=0; track<ntracks;track++) { // loop over tracks | |
197 | gAlice->ResetHits(); | |
198 | Int_t nbytes += TH->GetEvent(track); | |
199 | if (MUON) { | |
200 | for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); | |
201 | mHit; | |
202 | mHit=(AliMUONHit*)MUON->NextHit()) // loop over GEANT muon hits | |
203 | { | |
204 | Int_t nch = mHit->Chamber(); // chamber number | |
205 | ||
b6340d07 | 206 | |
207 | Int_t ipart = mHit->Particle(); | |
208 | Double_t x = mHit->X(); // x-pos of hit in chamber #1 | |
209 | Double_t y = mHit->Y(); // y-pos of hit in chamber #1 | |
210 | Double_t z = mHit->Z(); // z-pos of hit in chamber #1 | |
b45fd22b | 211 | |
212 | Int_t iSt = (nch-1)/2; | |
213 | if ((nch-1) < 10) { | |
214 | if (ipart == kMuonPlus){ | |
215 | hitChamberMUPlus[nch-1] = kTRUE; | |
216 | hitStationMUPlus[iSt] = kTRUE; | |
217 | } | |
218 | if (ipart == kMuonMinus){ | |
219 | hitChamberMUMinus[nch-1] = kTRUE; | |
220 | hitStationMUMinus[iSt] = kTRUE; | |
221 | } | |
222 | nHitChamber++; | |
223 | } | |
b6340d07 | 224 | |
b45fd22b | 225 | if (nch != 1) continue; |
b6340d07 | 226 | |
227 | if (ipart == kMuonPlus || ipart == kMuonMinus) { | |
228 | Double_t px0=mHit->Px(); | |
229 | Double_t py0=mHit->Py(); | |
230 | Double_t pz0=mHit->Pz(); | |
231 | Double_t nonBendingSlope=px0/pz0; | |
232 | Double_t bendingSlope=py0/pz0; | |
233 | ||
234 | // create an AliMUONTrackParam object in chamber #1 | |
235 | AliMUONTrackParam *trackParam = new AliMUONTrackParam(); | |
b45fd22b | 236 | Int_t sign = -TMath::Sign(1.0,ipart); |
237 | Double_t bendingMometum = sign* TMath::Sqrt(pz0*pz0+py0*py0); // Bug px -> py !!! | |
b6340d07 | 238 | Double_t inverseBendingMomentum = 1/bendingMometum; |
239 | ||
240 | trackParam->SetBendingCoor(y); | |
241 | trackParam->SetNonBendingCoor(x); | |
242 | trackParam->SetInverseBendingMomentum(inverseBendingMomentum); | |
243 | trackParam->SetBendingSlope(bendingSlope); | |
244 | trackParam->SetNonBendingSlope(nonBendingSlope); | |
245 | trackParam->SetZ(z); | |
246 | ||
247 | if (printLevel >= 1) { | |
248 | cout <<" **** before extrap " <<endl; | |
249 | cout << " px, py, pz = " <<px0<<" "<<py0<<" "<<pz0<<endl; | |
250 | trackParam->Dump(); | |
251 | } | |
b45fd22b | 252 | if (!testSingle) |
b6340d07 | 253 | trackParam->ExtrapToVertex(); // extrapolate track parameters through the absorber |
254 | if (printLevel >= 1) { | |
255 | cout <<" **** after extrap " <<endl; | |
256 | trackParam->Dump(); | |
257 | } | |
258 | ||
259 | Int_t iMuon; | |
260 | if (ipart == kMuonPlus) iMuon = 0; | |
261 | else iMuon = 1; | |
262 | ||
263 | // calculate track parameters extrapolated to vertex. | |
264 | bendingSlope = trackParam->GetBendingSlope(); | |
265 | nonBendingSlope = trackParam->GetNonBendingSlope(); | |
266 | Double_t pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum()); | |
267 | pZ[iMuon] = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope); | |
268 | pX[iMuon] = pZ[iMuon] * nonBendingSlope; | |
269 | pY[iMuon] = pZ[iMuon] * bendingSlope; | |
270 | pTot[iMuon] = TMath::Sqrt(pYZ*pYZ+pX[iMuon]*pX[iMuon]); | |
b45fd22b | 271 | theta[iMuon] = TMath::ATan2(TMath::Sqrt(pX[iMuon]*pX[iMuon]+pY[iMuon]*pY[iMuon]),pZ[iMuon]); |
b6340d07 | 272 | radius[iMuon] = TMath::Sqrt(x*x+y*y); |
273 | ||
274 | if (printLevel >= 1) | |
275 | cout << " px, py, pz, theta, radius = " <<pX[iMuon]<<" "<<pY[iMuon]<<" "<<pZ[iMuon]<<" "<<theta[iMuon]<<" "<<radius[iMuon]<<endl; | |
276 | delete trackParam; | |
277 | ||
278 | } // if MuonPlus or MuonMinus | |
279 | } // loop over MUON hits | |
280 | } // if MUON | |
b45fd22b | 281 | } // loop over tracks |
282 | ||
283 | ||
284 | Bool_t goodMUPlus = kTRUE; | |
285 | Bool_t goodMUMinus = kTRUE; | |
286 | for (Int_t iCh=0; iCh < 1; iCh++) { // !!!! 10 -> 1 | |
287 | if (!hitChamberMUPlus[iCh]) { | |
288 | goodMUPlus = kFALSE; | |
289 | printf(" evt number %d no muon+ in chamber %d \n",nev,iCh+1); | |
290 | } | |
291 | if (!hitChamberMUMinus[iCh]) { | |
292 | goodMUMinus = kFALSE; | |
293 | printf(" evt number %d no muon- in chamber %d \n",nev,iCh+1); | |
294 | } | |
295 | } | |
296 | ||
297 | // for (Int_t iSt=0; iSt < 5; iSt++) { | |
298 | // if (!hitStationMUPlus[iSt]) { | |
299 | // goodMUPlus = kFALSE; | |
300 | // // printf(" evt number %d no muon+ in chamber %d \n",nev,iCh+1); | |
301 | // } | |
302 | // if (!hitStationMUMinus[iSt]) { | |
303 | // goodMUMinus = kFALSE; | |
304 | // // printf(" evt number %d no muon- in chamber %d \n",nev,iCh+1); | |
305 | // } | |
306 | // } | |
b6340d07 | 307 | |
308 | if (pX[0] != 0 && pX[1] != 0) { // if track 1 & 2 exist | |
309 | Float_t massResonance = MuPlusMuMinusMass(pX[0],pY[0],pZ[0],pX[1],pY[1],pZ[1],massMuon); | |
310 | hInvMassRes->Fill(massResonance); | |
b45fd22b | 311 | } |
312 | ||
313 | ||
314 | for (Int_t i = 0; i<2; i++){ // !!! 1 -> 2 | |
315 | if (i == 0 && !goodMUPlus) continue; | |
316 | if (i == 1 && !goodMUMinus) continue; | |
317 | if (pX[i] == 0) continue; | |
318 | ||
319 | hRap->Fill(RapGen[i]); | |
320 | ||
321 | Double_t cosAngle = pX[i]*PxGen[i]+pY[i]*PyGen[i]+pZ[i]*PzGen[i]; | |
322 | cosAngle = cosAngle/(pTot[i]*PTotGen[i]); | |
323 | Double_t deltaAngle = TMath::ACos(cosAngle)*180/TMath::Pi(); | |
324 | Double_t deltaPNorm = 0; | |
325 | if (testSingle) { | |
326 | deltaPNorm = (PTotGen[i]-pTot[i])* TMath::Cos(ThetaGen[i]); | |
327 | } else { | |
328 | deltaPNorm = PTotGen[i]-pTot[i]; | |
329 | } | |
330 | ||
331 | Float_t thetaGX = TMath::ATan2(PxGen[i],PzGen[i]); | |
332 | Float_t thetaRX = TMath::ATan2(pX[i],pZ[i]); | |
333 | Float_t deltaAngleX = thetaRX - thetaGX; | |
334 | ||
335 | Float_t thetaGY = TMath::ATan2(PyGen[i],PzGen[i]); | |
336 | Float_t thetaRY = TMath::ATan2(pY[i],pZ[i]); | |
337 | Float_t deltaAngleY = thetaRY - thetaGY; | |
338 | ||
339 | Float_t phiG = TMath::ATan2(PyGen[i],PxGen[i])*180/TMath::Pi(); | |
340 | Float_t phiR = TMath::ATan2(pY[i],pX[i])*180/TMath::Pi(); | |
341 | ||
342 | hDeltaPxvsPhi->Fill(phiG,PxGen[i]-pX[i]); | |
343 | hDeltaPyvsPhi->Fill(phiG,PyGen[i]-pY[i]); | |
344 | hDeltaPhivsPz->Fill(PzGen[i],phiR-phiG); | |
345 | ||
346 | Int_t iHist=0; | |
347 | if (ThetaGen[i] < (3*TMath::Pi()/180)) { | |
348 | hDeltaP1->Fill(pTot[i]-PTotGen[i]); | |
349 | hDeltaAngle1->Fill(deltaAngle); | |
350 | hP1->Fill(PTotGen[i]); | |
351 | hProfDeltaPvsP1->Fill(PTotGen[i],deltaPNorm); | |
352 | h2DeltaPvsP1->Fill(PTotGen[i],deltaPNorm); | |
353 | iHist = (PTotGen[i]-xProfMin)*nBinProf/(xProfMax-xProfMin); | |
354 | if (PTotGen[i] > xProfMin && PTotGen[i] < xProfMax ) { | |
355 | DeltaP1[iHist]->Fill(deltaPNorm); | |
356 | DeltaAngleX1[iHist]->Fill(deltaAngleX); | |
357 | DeltaAngleY1[iHist]->Fill(deltaAngleY); | |
358 | } | |
359 | } else { | |
360 | hDeltaP2->Fill(pTot[i]-PTotGen[i]); | |
361 | hDeltaAngle2->Fill(deltaAngle); | |
362 | hP2->Fill(PTotGen[i]); | |
363 | hProfDeltaPvsP2->Fill(PTotGen[i],deltaPNorm); | |
364 | h2DeltaPvsP2->Fill(PTotGen[i],deltaPNorm); | |
365 | iHist = (PTotGen[i]-xProfMin)*nBinProf/(xProfMax-xProfMin); | |
366 | if (PTotGen[i] > xProfMin && PTotGen[i] < xProfMax ) { | |
367 | DeltaP2[iHist]->Fill(deltaPNorm); | |
368 | DeltaAngleX2[iHist]->Fill(deltaAngleX); | |
369 | DeltaAngleY2[iHist]->Fill(deltaAngleY); | |
b6340d07 | 370 | } |
371 | } | |
372 | } | |
373 | } // loop over event | |
b45fd22b | 374 | |
375 | ||
376 | Float_t sigmaP = 0; | |
377 | Float_t xBin = 0; | |
378 | for (Int_t iBin = 1; iBin <= nBinProf; iBin++){ | |
379 | sigmaP = hProfDeltaPvsP1->GetBinError(iBin); | |
380 | xBin = hProfDeltaPvsP1->GetBinCenter(iBin); | |
381 | hSigmaPvsP1->SetBinContent(iBin,sigmaP/xBin); | |
382 | sigmaP = hProfDeltaPvsP2->GetBinError(iBin); | |
383 | xBin = hProfDeltaPvsP2->GetBinCenter(iBin); | |
384 | hSigmaPvsP2->SetBinContent(iBin,sigmaP/xBin); | |
385 | } | |
386 | ||
387 | TF1 *g1= new TF1("g1","gaus",-25,25) ; | |
388 | TF1 *g2= new TF1("g2","gaus",-25,25) ; | |
389 | Float_t sigmaPGaus; | |
390 | Float_t xBinGaus; | |
391 | for (Int_t iHist=0; iHist < nBinProf; iHist++){ | |
392 | DeltaP1[iHist]->Fit("g1","RQ"); | |
393 | sigmaPGaus = g1->GetParameter(2); | |
394 | printf(" ** 1 ** iHist= %d , sigmaPGaus = %f \n",iHist,sigmaPGaus); | |
395 | xBinGaus = hSigGausPvsP1->GetBinCenter(iHist+1); | |
396 | hSigGausPvsP1->SetBinContent(iHist+1,sigmaPGaus/xBinGaus); | |
397 | ||
398 | DeltaP2[iHist]->Fit("g2","RQ"); | |
399 | sigmaPGaus = g2->GetParameter(2); | |
400 | printf(" ** 2 ** iHist= %d , sigmaGaus = %f \n",iHist,sigmaPGaus); | |
401 | xBinGaus = hSigGausPvsP2->GetBinCenter(iHist+1); | |
402 | hSigGausPvsP2->SetBinContent(iHist+1,sigmaPGaus/xBinGaus); | |
403 | } | |
b6340d07 | 404 | |
b45fd22b | 405 | |
406 | // Angles | |
407 | TF1 *g3= new TF1("g3","gaus") ; | |
408 | TF1 *g4= new TF1("g4","gaus") ; | |
409 | TF1 *g5= new TF1("g5","gaus") ; | |
410 | TF1 *g6= new TF1("g6","gaus") ; | |
411 | Float_t sigmaAngleGaus,limitGaus,errorSigma; | |
412 | Float_t nSigFit = 3; | |
413 | for (Int_t iHist=0; iHist < nBinProf; iHist++){ | |
414 | limitGaus = nSigFit * (DeltaAngleX1[iHist]->GetRMS()); | |
415 | g3->SetRange(-limitGaus,limitGaus); | |
416 | DeltaAngleX1[iHist]->Fit("g3","RQ"); | |
417 | sigmaAngleGaus = g3->GetParameter(2); | |
418 | errorSigma = g3->GetParError(2); | |
419 | printf(" ** 1 ** iHist= %d , sigmaAngleGaus X = %f \n",iHist,sigmaAngleGaus); | |
420 | hSigGausAngleXvsP1->SetBinContent(iHist+1,sigmaAngleGaus); | |
421 | hSigGausAngleXvsP1->SetBinError(iHist+1,errorSigma); | |
422 | ||
423 | limitGaus = nSigFit * (DeltaAngleY1[iHist]->GetRMS()); | |
424 | g4->SetRange(-limitGaus,limitGaus); | |
425 | DeltaAngleY1[iHist]->Fit("g4","RQ"); | |
426 | sigmaAngleGaus = g4->GetParameter(2); | |
427 | errorSigma = g4->GetParError(2); | |
428 | printf(" ** 1 ** iHist= %d , sigmaAngleGaus Y = %f \n",iHist,sigmaAngleGaus); | |
429 | hSigGausAngleYvsP1->SetBinContent(iHist+1,sigmaAngleGaus); | |
430 | hSigGausAngleYvsP1->SetBinError(iHist+1,errorSigma); | |
431 | ||
432 | limitGaus = nSigFit * (DeltaAngleX2[iHist]->GetRMS()); | |
433 | g5->SetRange(-limitGaus,limitGaus); | |
434 | DeltaAngleX2[iHist]->Fit("g5","RQ"); | |
435 | sigmaAngleGaus = g5->GetParameter(2); | |
436 | errorSigma = g5->GetParError(2); | |
437 | printf(" ** 1 ** iHist= %d , sigmaAngleGaus X = %f \n",iHist,sigmaAngleGaus); | |
438 | hSigGausAngleXvsP2->SetBinContent(iHist+1,sigmaAngleGaus); | |
439 | hSigGausAngleXvsP2->SetBinError(iHist+1,errorSigma); | |
440 | ||
441 | limitGaus = nSigFit * (DeltaAngleY2[iHist]->GetRMS()); | |
442 | g6->SetRange(-limitGaus,limitGaus); | |
443 | DeltaAngleY2[iHist]->Fit("g6","RQ"); | |
444 | sigmaAngleGaus = g6->GetParameter(2); | |
445 | errorSigma = g6->GetParError(2); | |
446 | printf(" ** 1 ** iHist= %d , sigmaAngleGaus Y = %f \n",iHist,sigmaAngleGaus); | |
447 | hSigGausAngleYvsP2->SetBinContent(iHist+1,sigmaAngleGaus); | |
448 | hSigGausAngleYvsP2->SetBinError(iHist+1,errorSigma); | |
449 | ||
450 | } | |
b6340d07 | 451 | // save histograms in MUONTestAbso.root |
452 | TFile *histoFile = new TFile("MUONTestAbso.root", "RECREATE"); | |
b45fd22b | 453 | |
b6340d07 | 454 | hInvMassRes->Write(); |
b45fd22b | 455 | hRap->Write(); |
456 | ||
b6340d07 | 457 | hDeltaP1->Write(); |
458 | hDeltaAngle1->Write(); | |
b45fd22b | 459 | hP1->Write(); |
460 | hProfDeltaPvsP1->Write(); | |
461 | h2DeltaPvsP1->Write(); | |
462 | hSigmaPvsP1->Write(); | |
463 | hSigGausPvsP1->Write(); | |
464 | hSigGausAngleXvsP1->Write(); | |
465 | hSigGausAngleYvsP1->Write(); | |
466 | ||
b6340d07 | 467 | hDeltaP2->Write(); |
468 | hDeltaAngle2->Write(); | |
b45fd22b | 469 | hP2->Write(); |
470 | hProfDeltaPvsP2->Write(); | |
471 | h2DeltaPvsP2->Write(); | |
472 | hSigmaPvsP2->Write(); | |
473 | hSigGausPvsP2->Write(); | |
474 | hSigGausAngleXvsP2->Write(); | |
475 | hSigGausAngleYvsP2->Write(); | |
476 | ||
477 | for (Int_t iHist=0; iHist < nBinProf; iHist++){ | |
478 | DeltaP1[iHist]->Write(); | |
479 | DeltaP2[iHist]->Write(); | |
480 | DeltaAngleX1[iHist]->Write(); | |
481 | DeltaAngleY1[iHist]->Write(); | |
482 | DeltaAngleX2[iHist]->Write(); | |
483 | DeltaAngleY2[iHist]->Write(); | |
484 | } | |
485 | ||
486 | hDeltaPxvsPhi->Write(); | |
487 | hDeltaPyvsPhi->Write(); | |
488 | hDeltaPhivsPz->Write(); | |
489 | ||
b6340d07 | 490 | histoFile->Close(); |
491 | ||
492 | } | |
493 | ||
494 | ||
495 | Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t MuonMass) | |
496 | { | |
497 | // return invariant mass for particle 1 & 2 whose masses are equal to MuonMass | |
498 | ||
499 | Float_t e1 = TMath::Sqrt(MuonMass * MuonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1); | |
500 | Float_t e2 = TMath::Sqrt(MuonMass * MuonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2); | |
501 | return (TMath::Sqrt(2.0 * (MuonMass * MuonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2))); | |
502 | } | |
b45fd22b | 503 | |
504 | Float_t rapParticle(Float_t Px, Float_t Py, Float_t Pz, Float_t Mass) | |
505 | { | |
506 | // return rapidity for particle | |
507 | ||
508 | // Particle energy | |
509 | Float_t Energy = TMath::Sqrt(Px*Px+Py*Py+Pz*Pz+Mass*Mass); | |
510 | // Rapidity | |
511 | Float_t Rapidity = 0.5*TMath::Log((Energy+Pz)/(Energy-Pz)); | |
512 | return Rapidity; | |
513 | } |