]>
Commit | Line | Data |
---|---|---|
bd022fc7 | 1 | #if !defined(__CINT__) || defined(__MAKECINT__) |
2 | ||
3 | // ROOT includes | |
4 | #include "TBranch.h" | |
5 | #include "TClonesArray.h" | |
6 | #include "TTree.h" | |
7 | #include "TNtuple.h" | |
8 | #include "TParticle.h" | |
9 | #include "TFile.h" | |
10 | ||
11 | #include "TH1.h" | |
12 | #include "TH1F.h" | |
13 | #include "TH2F.h" | |
14 | #include "TF1.h" | |
15 | #include "TMath.h" | |
16 | ||
17 | #include "TCanvas.h" | |
18 | #include "TGraph.h" | |
19 | #include "TGraphErrors.h" | |
20 | #include "TGraph2D.h" | |
21 | #include "TStyle.h" | |
22 | #include "TFitter.h" | |
23 | #include "TRandom.h" | |
24 | ||
25 | // STEER includes | |
26 | #include "AliRun.h" | |
27 | #include "AliRunLoader.h" | |
28 | #include "AliHeader.h" | |
29 | #include "AliLoader.h" | |
30 | #include "AliTracker.h" | |
31 | #include "AliStack.h" | |
32 | #include "AliMagFMaps.h" | |
33 | ||
34 | ||
35 | // MUON includes | |
36 | #include "AliMUON.h" | |
37 | #include "AliMUONData.h" | |
38 | #include "AliMUONConstants.h" | |
39 | ||
40 | #include "AliMUONHit.h" | |
41 | #include "AliMUONHitForRec.h" | |
42 | #include "AliMUONRawCluster.h" | |
43 | #include "AliMUONTrack.h" | |
44 | #include "AliMUONTrackParam.h" | |
45 | #include "AliMUONTrackExtrap.h" | |
46 | ||
47 | #include "AliMpVSegmentation.h" | |
48 | #include "AliMpIntPair.h" | |
49 | #include "AliMpDEManager.h" | |
50 | #endif | |
51 | ||
52 | ||
53 | //Macro to calculate the resolution and the efficiency of chamber(s) chosen in function | |
54 | //of Phi (angle on the chamber between the X axis and the straight line created by the | |
55 | //center of the chamber and the impact of particles on this chamber, the azimuthal angle) | |
56 | //and Theta (the polar angle), or in function of ThetaI (angle of incidence of particles | |
57 | //on the chamber) | |
58 | ||
59 | ||
60 | ||
61 | ||
62 | const Double_t kInvPi = 1./3.14159265; | |
63 | ||
64 | ||
65 | //Chamber number: | |
66 | Int_t chamberNbr; | |
67 | //Number of events: | |
68 | Int_t nEvents, iEvent; | |
69 | //Number of tracks: | |
70 | Int_t nTracks, iTrack; | |
71 | //Number of hits: | |
72 | Int_t nHits,iHit; | |
73 | //Chamber(s) chosen | |
74 | Int_t firstChamber, lastChamber; | |
75 | ||
76 | ||
77 | AliMUONTrack * track ; | |
78 | AliMUONHitForRec * hit = 0; | |
79 | AliMUONTrackParam * trackParam = 0; | |
80 | ||
81 | TClonesArray * tracks ; | |
82 | TClonesArray * trackParams; | |
83 | TClonesArray * hits ; | |
84 | ||
85 | ||
86 | ||
87 | ||
88 | ||
89 | ||
90 | /*****************************************************************************************************************/ | |
91 | /*****************************************************************************************************************/ | |
92 | /*EFFICIENCY*/ | |
93 | ||
94 | void efficiency( Int_t event2Check=0, char * filename="galice.root" ) | |
95 | { | |
96 | cout<<"\nChamber(s) chosen;\nFirst chamber:"; | |
97 | cin>>firstChamber; | |
98 | cout<<"Last chamber:"; | |
99 | cin>>lastChamber; | |
100 | cout<<"\n\n"; | |
101 | ||
102 | //Creating Run Loader and openning file containing Hits | |
103 | //-------------------------------------------------------------------------------------------------------------- | |
104 | AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ"); | |
105 | if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;} | |
106 | ||
107 | ||
108 | //Getting MUONLoader | |
109 | //-------------------------------------------------------------------------------------------------------------- | |
110 | AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader"); | |
111 | MUONLoader -> LoadTracks("READ"); | |
112 | MUONLoader -> LoadRecPoints("READ"); | |
113 | ||
114 | ||
115 | //Creating a MUON data container | |
116 | //-------------------------------------------------------------------------------------------------------------- | |
117 | AliMUONData muondata(MUONLoader,"MUON","MUON"); | |
118 | ||
119 | ||
120 | nEvents = (Int_t) RunLoader -> GetNumberOfEvents(); | |
121 | ||
122 | //Loop on events | |
123 | Int_t trackNb = 0; | |
124 | Int_t chamber[10] = {0}; | |
125 | Int_t detEltNew, detElt; | |
126 | Int_t detEltOld = 0; | |
127 | ||
128 | for ( iEvent = 0; iEvent < nEvents; ++iEvent ) | |
129 | { | |
130 | if ( event2Check!=0 ) | |
131 | iEvent=event2Check; | |
132 | printf("\r>>> Event %d ",iEvent); | |
133 | RunLoader -> GetEvent(iEvent); | |
134 | //Addressing | |
135 | muondata.SetTreeAddress("RT"); | |
136 | muondata.GetRecTracks(); | |
137 | tracks = muondata.RecTracks(); | |
138 | ||
139 | //Loop on track | |
140 | nTracks = (Int_t) tracks -> GetEntriesFast(); | |
141 | for ( iTrack = 0; iTrack < nTracks; ++iTrack ) | |
142 | { | |
143 | track = (AliMUONTrack*) tracks -> At(iTrack); | |
144 | hits = track -> GetHitForRecAtHit(); | |
145 | detEltOld = 0; | |
146 | //Loop on hit | |
147 | nHits = (Int_t) hits -> GetEntriesFast(); | |
148 | ||
149 | for ( iHit = 0; iHit < nHits; ++iHit ) | |
150 | { | |
151 | hit = (AliMUONHitForRec*) hits -> At(iHit); | |
152 | chamberNbr = hit -> GetChamberNumber(); | |
153 | detElt = hit -> GetDetElemId(); | |
154 | detEltNew = int(detElt/100); | |
155 | if( chamberNbr >= firstChamber-1 ) { | |
156 | if( chamberNbr < lastChamber ) { | |
157 | if( detEltNew == detEltOld ) | |
158 | continue; | |
159 | else { | |
160 | chamber[chamberNbr] = chamber[chamberNbr] + 1; | |
161 | detEltOld = detEltNew; | |
162 | } | |
163 | } | |
164 | } | |
165 | } | |
166 | //End loop on hit | |
167 | ||
168 | } | |
169 | //End loop on track | |
170 | muondata.ResetRecTracks(); | |
171 | if (event2Check != 0) | |
172 | iEvent = nEvents; | |
173 | trackNb = trackNb + nTracks; | |
174 | } | |
175 | //End loop on event | |
176 | //-------------------------------------------------------------------------------------------------------------- | |
177 | ||
178 | cout<<"\n\n\n"; | |
179 | for (Int_t i = firstChamber-1; i < lastChamber; ++i ) | |
180 | { | |
181 | printf ("\nChamber %d has responded %d times on %d tracks", i+1, chamber[i], trackNb); | |
182 | if (trackNb == 0) | |
183 | cout<<"\nEfficiency = ? (IS UNKNOWN)\n"; | |
184 | else { | |
185 | Double_t eff = chamber[i]*100./trackNb; cout<<"\nEfficiency = "<<eff<<" %\n"; | |
186 | } | |
187 | } | |
188 | cout<<"\n\n\n"; | |
189 | MUONLoader->UnloadTracks(); | |
190 | } | |
191 | ||
192 | ||
193 | ||
194 | ||
195 | ||
196 | /*****************************************************************************************************************/ | |
197 | /*****************************************************************************************************************/ | |
198 | /*EFFICIENCY IN FUNCTION OF THETA AND PHI*/ | |
199 | ||
200 | void efficiencyThetaPhi( Int_t event2Check=0, char * filename="galice.root" ) | |
201 | { | |
202 | cout<<"\nChamber(s) chosen;\nFirst chamber:"; | |
203 | cin>>firstChamber; | |
204 | cout<<"Last chamber:"; | |
205 | cin>>lastChamber; | |
206 | cout<<"\n\n"; | |
207 | ||
208 | Int_t eff [54] = {0}; | |
209 | Int_t trackN [54] = {0}; | |
210 | Int_t chamber; | |
211 | Int_t oldChamber; | |
212 | ||
213 | //Creating Run Loader and openning file containing Hits | |
214 | //-------------------------------------------------------------------------------------------------------------- | |
215 | AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ"); | |
216 | if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;} | |
217 | ||
218 | ||
219 | //Getting MUONLoader | |
220 | //-------------------------------------------------------------------------------------------------------------- | |
221 | AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader"); | |
222 | MUONLoader -> LoadTracks("READ"); | |
223 | MUONLoader -> LoadRecPoints("READ"); | |
224 | ||
225 | ||
226 | //-------------------------------------------------------------------------------------------------------------- | |
227 | //Set mag field; waiting for mag field in CDB | |
228 | printf("Loading field map...\n"); | |
229 | if (!AliTracker::GetFieldMap()) { | |
230 | AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG); | |
231 | AliTracker::SetFieldMap(field, kFALSE);} | |
232 | ||
233 | ||
234 | //-------------------------------------------------------------------------------------------------------------- | |
235 | //Set Field Map for track extrapolation | |
236 | AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap()); | |
237 | ||
238 | ||
239 | //Creating a MUON data container | |
240 | //-------------------------------------------------------------------------------------------------------------- | |
241 | AliMUONData muondata(MUONLoader,"MUON","MUON"); | |
242 | ||
243 | ||
244 | nEvents = (Int_t) RunLoader -> GetNumberOfEvents(); | |
245 | ||
246 | //Loop on events | |
247 | for ( iEvent = 0; iEvent < nEvents; ++iEvent ) | |
248 | { | |
249 | if ( event2Check!=0 ) | |
250 | iEvent=event2Check; | |
251 | printf("\r>>> Event %d ",iEvent); | |
252 | RunLoader->GetEvent(iEvent); | |
253 | ||
254 | //Addressing | |
255 | muondata.SetTreeAddress("RT"); | |
256 | muondata.GetRecTracks(); | |
257 | tracks = muondata.RecTracks(); | |
258 | ||
259 | //Loop on track | |
260 | nTracks = (Int_t) tracks -> GetEntriesFast(); | |
261 | for ( iTrack = 0; iTrack < nTracks; ++iTrack ) | |
262 | { | |
263 | track = (AliMUONTrack*) tracks ->At(iTrack) ; | |
264 | trackParams = track ->GetTrackParamAtHit(); | |
265 | hits = track ->GetHitForRecAtHit() ; | |
266 | chamber = firstChamber-1; | |
267 | oldChamber = -1; | |
268 | Int_t k = 1; | |
269 | ||
270 | //Loop on hits | |
271 | nHits = (Int_t) hits->GetEntriesFast(); | |
272 | for ( iHit = 0; iHit<nHits; ++iHit ) | |
273 | { | |
274 | trackParam = (AliMUONTrackParam*) trackParams -> At(iHit); | |
275 | hit = (AliMUONHitForRec* ) hits -> At(iHit); | |
276 | chamberNbr = hit -> GetChamberNumber(); | |
277 | ||
278 | if ( chamberNbr == oldChamber ) | |
279 | continue; | |
280 | else { | |
281 | oldChamber = chamberNbr; | |
282 | if ( chamberNbr > chamber - k ) { | |
283 | if ( chamber < lastChamber ) { | |
284 | if ( chamberNbr == chamber ) { | |
285 | //Positions | |
286 | Double_t traX, traY, traZ; | |
287 | Double_t hitX, hitY, hitZ; | |
288 | Double_t aveX, aveY ; | |
289 | ||
290 | //Angle (Phi) | |
291 | Double_t phi = 0.; | |
292 | Double_t theta = 0.; | |
293 | Int_t iPhi = 0 ; | |
294 | Int_t iTheta = 0 ; | |
295 | ||
296 | traX = trackParam -> GetNonBendingCoor(); | |
297 | traY = trackParam -> GetBendingCoor() ; | |
298 | traZ = trackParam -> GetZ() ; | |
299 | ||
300 | hitX = hit -> GetNonBendingCoor(); | |
301 | hitY = hit -> GetBendingCoor() ; | |
302 | hitZ = hit -> GetZ() ; | |
303 | ||
304 | aveX = 1./2.*(traX + hitX); | |
305 | aveY = 1./2.*(traY + hitY); | |
306 | ||
307 | //The calculation of phi: | |
308 | phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX )); | |
309 | ||
310 | //The calculation of theta, theta is in fact 180 - Theta(The polar angle): | |
311 | theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ )); | |
312 | ||
313 | if ( phi < 0.) phi = 360 - abs(phi); | |
314 | iPhi = int( phi/72. ); | |
315 | iTheta = int( theta ); | |
316 | if( theta > 10 ) iTheta = 9; | |
317 | if( theta < 1 ) iTheta = 1; | |
318 | ||
319 | eff [9*iPhi+iTheta-1] = eff [9*iPhi+iTheta-1] + 1; | |
320 | trackN [9*iPhi+iTheta-1] = trackN [9*iPhi+iTheta-1] + 1; | |
321 | chamber = chamber + 1; | |
322 | } | |
323 | ||
324 | else { | |
325 | //Positions | |
326 | Double_t chamberZpos; | |
327 | Double_t exXpos, exYpos; | |
328 | ||
329 | //Angles | |
330 | Double_t phi = 0.; | |
331 | Double_t theta = 0.; | |
332 | Int_t iPhi = 0 ; | |
333 | Int_t iTheta = 0 ; | |
334 | ||
335 | chamberZpos = AliMUONConstants::DefaultChamberZ(chamber); | |
336 | AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos); | |
337 | exXpos = (Double_t) trackParam->GetNonBendingCoor(); | |
338 | exYpos = (Double_t) trackParam->GetBendingCoor(); | |
339 | ||
340 | //The calculation of phi: | |
341 | phi = 180. * kInvPi * (TMath::ATan2( exYpos, exXpos )); | |
342 | ||
343 | //The calculation of theta, theta is in fact 180 - Theta(The polar angle): | |
344 | theta = 180. * kInvPi * (TMath::ATan2( sqrt(exXpos*exXpos+exYpos*exYpos), - chamberZpos )); | |
345 | ||
346 | if ( phi < 0.) phi = 360 - abs(phi); | |
347 | iPhi = int( phi/72. ); | |
348 | iTheta = int( theta ); | |
349 | if( theta > 10 ) iTheta = 9; | |
350 | if( theta < 1 ) iTheta = 1; | |
351 | ||
352 | eff [9*iPhi+iTheta-1] = eff [9*iPhi+iTheta-1] + 0; | |
353 | trackN [9*iPhi+iTheta-1] = trackN [9*iPhi+iTheta-1] + 1; | |
354 | chamber = chamber + 1; | |
355 | iHit = iHit - 1; | |
356 | } | |
357 | } | |
358 | } | |
359 | } | |
360 | ||
361 | if ( iHit == nHits-1 ) { | |
362 | if ( chamber == lastChamber ) | |
363 | continue; | |
364 | else { | |
365 | oldChamber = -1; | |
366 | k = 5; | |
367 | iHit = iHit-1; | |
368 | } | |
369 | } | |
370 | } | |
371 | //End Loop on hits | |
372 | ||
373 | } | |
374 | //End Loop on tracks | |
375 | ||
376 | muondata.ResetRecTracks(); | |
377 | if ( event2Check!=0 ) | |
378 | iEvent=nEvents; | |
379 | } | |
380 | //End Loop on events | |
381 | ||
382 | ||
383 | TCanvas * c1 = new TCanvas(); | |
384 | TGraph2D* effPhiTheta = new TGraph2D(); | |
385 | Double_t efficiency = 0; | |
386 | ||
387 | if ( firstChamber == lastChamber ) | |
388 | { | |
389 | for ( Int_t ph = 0; ph < 5; ++ph ) | |
390 | { | |
391 | for ( Int_t th = 1; th < 10; ++th ) | |
392 | { | |
393 | Int_t i = 9*ph+th-1; | |
394 | cout<<"\nFor Phi = ["<<ph*72<<","<<ph*72+72<<"] and Theta = ["<<179-th<<","<<180-th<<"]:"; | |
395 | cout<<"\nThe chamber "<<firstChamber<<" has responded "<<eff[i]<<" times on "<<trackN[i]<<" tracks\n"; | |
396 | ||
397 | Double_t e = eff [i] ; | |
398 | Double_t n = trackN[i] ; | |
399 | Double_t p = ph*72.+36.; | |
400 | Double_t t = th*1. +0.5; | |
401 | ||
402 | if ( trackN[i] == 0 ) { | |
403 | efficiency = 0.; | |
404 | cout<<"Efficiency = ? % ( IS UNKNOWN )\n"; | |
405 | } | |
406 | else { | |
407 | efficiency = 100.*e/n; | |
408 | cout<<"Efficiency = "<<efficiency<<" %\n"; | |
409 | } | |
410 | ||
411 | effPhiTheta -> SetPoint( i, p, t, efficiency); | |
412 | } | |
413 | } | |
414 | } | |
415 | ||
416 | else{ | |
417 | for ( Int_t ph = 0; ph < 5; ++ph ) | |
418 | { | |
419 | for ( Int_t th = 1; th < 10; ++th ) | |
420 | { | |
421 | Int_t i = 9*ph+th-1; | |
422 | cout<<"\nFor Phi = ["<<ph*72<<","<<ph*72+72<<"] and Theta = ["<<179-th<<","<<180-th<<"]:"; | |
423 | cout<<"\nThe chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff[i]<<" times on "<<trackN[i]<<" tracks\n"; | |
424 | ||
425 | Double_t e = eff [i] ; | |
426 | Double_t n = trackN[i] ; | |
427 | Double_t p = ph*72.+36.; | |
428 | Double_t t = th*1. +0.5; | |
429 | ||
430 | if ( trackN[i] == 0 ) { | |
431 | efficiency = 0.; | |
432 | cout<<"Efficiency = ? % ( IS UNKNOWN )\n"; | |
433 | } | |
434 | else { | |
435 | efficiency = 100.*e/n; | |
436 | cout<<"Efficiency = "<<efficiency<<" %\n"; | |
437 | } | |
438 | ||
439 | effPhiTheta -> SetPoint( i, p, t, efficiency); | |
440 | } | |
441 | } | |
442 | } | |
443 | ||
444 | gStyle->SetPalette(1); | |
445 | effPhiTheta -> Draw("surf1"); | |
446 | ||
447 | cout<<"\n\n\n"; | |
448 | MUONLoader->UnloadTracks(); | |
449 | c1->Update(); | |
450 | ||
451 | } | |
452 | ||
453 | ||
454 | ||
455 | ||
456 | ||
457 | ||
458 | ||
459 | /*****************************************************************************************************************/ | |
460 | /*****************************************************************************************************************/ | |
461 | /*EFFICIENCY IN FUNCTION OF THETA OF INCIDENCE*/ | |
462 | ||
463 | ||
464 | void efficiencyThetaI( Int_t event2Check=0, char * filename="galice.root" ) | |
465 | { | |
466 | cout<<"\nChamber(s) chosen;\nFirst chamber:"; | |
467 | cin>>firstChamber; | |
468 | cout<<"Last chamber:"; | |
469 | cin>>lastChamber; | |
470 | cout<<"\n\n"; | |
471 | ||
472 | Int_t eff [12] = {0}; | |
473 | Int_t trackN [12] = {0}; | |
474 | Int_t chamber; | |
475 | Int_t oldChamber; | |
476 | ||
477 | ||
478 | //Creating Run Loader and openning file containing Hits | |
479 | //-------------------------------------------------------------------------------------------------------------- | |
480 | AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ"); | |
481 | if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;} | |
482 | ||
483 | ||
484 | //Getting MUONLoader | |
485 | //-------------------------------------------------------------------------------------------------------------- | |
486 | AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader"); | |
487 | MUONLoader -> LoadTracks("READ"); | |
488 | MUONLoader -> LoadRecPoints("READ"); | |
489 | ||
490 | ||
491 | //-------------------------------------------------------------------------------------------------------------- | |
492 | //Set mag field; waiting for mag field in CDB | |
493 | printf("Loading field map...\n"); | |
494 | if (!AliTracker::GetFieldMap()) { | |
495 | AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG); | |
496 | AliTracker::SetFieldMap(field, kFALSE);} | |
497 | ||
498 | ||
499 | //-------------------------------------------------------------------------------------------------------------- | |
500 | //Set Field Map for track extrapolation | |
501 | AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap()); | |
502 | ||
503 | ||
504 | //Creating a MUON data container | |
505 | //-------------------------------------------------------------------------------------------------------------- | |
506 | AliMUONData muondata(MUONLoader,"MUON","MUON"); | |
507 | ||
508 | ||
509 | nEvents = (Int_t) RunLoader -> GetNumberOfEvents(); | |
510 | ||
511 | //Loop on events | |
512 | for ( iEvent = 0; iEvent < nEvents; ++iEvent ) | |
513 | { | |
514 | if ( event2Check!=0 ) | |
515 | iEvent=event2Check; | |
516 | printf("\r>>> Event %d ",iEvent); | |
517 | RunLoader->GetEvent(iEvent); | |
518 | ||
519 | //Addressing | |
520 | muondata.SetTreeAddress("RT"); | |
521 | muondata.GetRecTracks(); | |
522 | tracks = muondata.RecTracks(); | |
523 | ||
524 | //Loop on track | |
525 | nTracks = (Int_t) tracks -> GetEntriesFast(); | |
526 | for ( iTrack = 0; iTrack < nTracks; ++iTrack ) | |
527 | { | |
528 | track = (AliMUONTrack*) tracks ->At(iTrack) ; | |
529 | trackParams = track ->GetTrackParamAtHit(); | |
530 | hits = track ->GetHitForRecAtHit() ; | |
531 | chamber = firstChamber - 1; | |
532 | oldChamber = -1; | |
533 | Int_t k = 1; | |
534 | ||
535 | //Loop on hits | |
536 | nHits = (Int_t) hits -> GetEntriesFast(); | |
537 | for ( iHit = 0; iHit < nHits; ++iHit ) | |
538 | { | |
539 | trackParam = (AliMUONTrackParam*) trackParams -> At(iHit); | |
540 | hit = (AliMUONHitForRec*) hits -> At(iHit); | |
541 | chamberNbr = hit -> GetChamberNumber(); | |
542 | ||
543 | if ( chamberNbr == oldChamber ) | |
544 | continue; | |
545 | else { | |
546 | oldChamber = chamberNbr; | |
547 | if ( chamberNbr > chamber - k ) { | |
548 | if ( chamber < lastChamber ) { | |
549 | if ( chamberNbr == chamber ) { | |
550 | //Momentum | |
551 | Double_t Px, Py, Pz, Pr; | |
552 | ||
553 | //Angle | |
554 | Double_t theta; | |
555 | Int_t iTheta; | |
556 | ||
557 | Px = trackParam -> Px() ; | |
558 | Py = trackParam -> Py() ; | |
559 | Pz = trackParam -> Pz() ; | |
560 | Pr = sqrt( Px*Px + Py*Py ); | |
561 | ||
562 | //The calculation of theta, the angle of incidence: | |
563 | theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr)); | |
564 | ||
565 | if ( theta < 79 ) iTheta = 11; | |
566 | else if ( theta < 90 ) iTheta = int( theta - 79.); | |
567 | else iTheta = 11; | |
568 | ||
569 | eff [iTheta] = eff [iTheta] + 1; | |
570 | trackN [iTheta] = trackN [iTheta] + 1; | |
571 | chamber = chamber + 1; | |
572 | } | |
573 | ||
574 | else { | |
575 | //Positions | |
576 | Double_t chamberZpos; | |
577 | ||
578 | //Momentum | |
579 | Double_t Px, Py, Pz, Pr; | |
580 | ||
581 | //Angles | |
582 | Double_t theta = 0.; | |
583 | Int_t iTheta = 0 ; | |
584 | ||
585 | chamberZpos = AliMUONConstants::DefaultChamberZ(chamber); | |
586 | AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos); | |
587 | ||
588 | Px = trackParam -> Px() ; | |
589 | Py = trackParam -> Py() ; | |
590 | Pz = trackParam -> Pz() ; | |
591 | Pr = sqrt( Px*Px + Py*Py ); | |
592 | ||
593 | //The calculation of thetaI, the angle of incidence: | |
594 | theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr)); | |
595 | ||
596 | if ( theta < 79 ) iTheta = 11; | |
597 | else if ( theta < 90 ) iTheta = int( theta - 79.); | |
598 | else iTheta = 11; | |
599 | ||
600 | eff [iTheta] = eff [iTheta] + 0; | |
601 | trackN [iTheta] = trackN [iTheta] + 1; | |
602 | chamber = chamber + 1; | |
603 | iHit = iHit - 1; | |
604 | } | |
605 | } | |
606 | } | |
607 | } | |
608 | ||
609 | if ( iHit == nHits-1 ) { | |
610 | if ( chamber == lastChamber ) | |
611 | continue; | |
612 | else { | |
613 | oldChamber = -1; | |
614 | k = 5; | |
615 | iHit = iHit-1; | |
616 | } | |
617 | } | |
618 | } | |
619 | //End loop on hits | |
620 | ||
621 | } | |
622 | //End Loop on tracks | |
623 | ||
624 | muondata.ResetRecTracks(); | |
625 | if ( event2Check!=0 ) | |
626 | iEvent=nEvents; | |
627 | } | |
628 | //End Loop on events | |
629 | ||
630 | Double_t t [11]; | |
631 | Double_t efficiency [11]; | |
632 | Int_t i = 11; | |
633 | ||
634 | Int_t th; | |
635 | TGraph * effTheta = new TGraph (); | |
636 | ||
637 | if ( firstChamber == lastChamber ) { | |
638 | for ( th = 0; th < 11; ++th ) | |
639 | { | |
640 | cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n"; | |
641 | cout<<"The chamber "<<firstChamber<<" has responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n"; | |
642 | ||
643 | t [th] = th + 79.5 ; | |
644 | Double_t e = eff [th]; | |
645 | Double_t n = trackN [th]; | |
646 | ||
647 | if ( n == 0. ) { | |
648 | efficiency [th] = 0.; | |
649 | cout<<"Efficiency = ? % (IS UNKNOWN) \n"; | |
650 | } | |
651 | else { | |
652 | efficiency [th] = 100.*e/n; | |
653 | cout<<"Efficiency = "<<efficiency [th]<<" %\n"; | |
654 | } | |
655 | } | |
656 | } | |
657 | ||
658 | else{ | |
659 | for ( th = 0; th < 11; ++th ) | |
660 | { | |
661 | cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n"; | |
662 | cout<<"The chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n"; | |
663 | ||
664 | t [th] = th + 79.5 ; | |
665 | Double_t e = eff [th]; | |
666 | Double_t n = trackN [th]; | |
667 | ||
668 | if ( n == 0. ) { | |
669 | efficiency [th] = 0.; | |
670 | cout<<"Efficiency = ? % (IS UNKNOWN) \n"; | |
671 | } | |
672 | else { | |
673 | efficiency [th] = 100.*e/n; | |
674 | cout<<"Efficiency = "<<efficiency [th]<<" %\n"; | |
675 | } | |
676 | } | |
677 | } | |
678 | ||
679 | TCanvas * c1 = new TCanvas (); | |
680 | effTheta = new TGraph( i, t, efficiency ); | |
681 | ||
682 | effTheta -> Draw("ALP"); | |
683 | ||
684 | MUONLoader->UnloadTracks(); | |
685 | c1->Update(); | |
686 | } | |
687 | ||
688 | ||
689 | ||
690 | ||
691 | ||
692 | ||
693 | /*****************************************************************************************************************/ | |
694 | /*****************************************************************************************************************/ | |
695 | /*RESOLUTION*/ | |
696 | ||
697 | void resolution( Int_t event2Check=0, char * filename="galice.root" ) | |
698 | { | |
699 | cout<<"\nChamber(s) chosen;\nFirst chamber:"; | |
700 | cin>>firstChamber; | |
701 | cout<<"Last chamber:"; | |
702 | cin>>lastChamber; | |
703 | cout<<"\n\n"; | |
704 | ||
705 | TH1F * hDeltax; | |
706 | TH1F * hDeltay; | |
707 | TH2 * hDelta3D; | |
708 | ||
709 | hDeltax = new TH1F ("hDeltax " , "No Interest", 100, -0.3, 0.3); | |
710 | hDeltay = new TH1F ("hDeltay " , "No Interest", 500, -0.1, 0.1); | |
711 | hDelta3D = new TH2F ("hDelta3D" , "No Interest", 100, -0.3, 0.3, 500, -0.1, 0.1); | |
712 | ||
713 | ||
714 | //Creating Run Loader and openning file containing Hits | |
715 | //-------------------------------------------------------------------------------------------------------------- | |
716 | AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ"); | |
717 | if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;} | |
718 | ||
719 | ||
720 | //Getting MUONLoader | |
721 | //-------------------------------------------------------------------------------------------------------------- | |
722 | AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader"); | |
723 | MUONLoader -> LoadTracks("READ"); | |
724 | MUONLoader -> LoadRecPoints("READ"); | |
725 | ||
726 | ||
727 | //Creating a MUON data container | |
728 | //-------------------------------------------------------------------------------------------------------------- | |
729 | AliMUONData muondata(MUONLoader,"MUON","MUON"); | |
730 | ||
731 | ||
732 | nEvents = (Int_t) RunLoader -> GetNumberOfEvents(); | |
733 | ||
734 | //Loop on events | |
735 | for ( iEvent = 0; iEvent < nEvents; ++iEvent ) | |
736 | { | |
737 | if (event2Check!=0) | |
738 | iEvent=event2Check; | |
739 | printf("\r>>> Event %d ",iEvent); | |
740 | RunLoader->GetEvent(iEvent); | |
741 | ||
742 | //Addressing | |
743 | muondata.SetTreeAddress("RT"); | |
744 | muondata.GetRecTracks(); | |
745 | tracks = muondata.RecTracks(); | |
746 | ||
747 | //Loop on track | |
748 | nTracks = (Int_t) tracks -> GetEntriesFast(); | |
749 | for ( iTrack = 0; iTrack < nTracks; ++iTrack ) | |
750 | { | |
751 | track = (AliMUONTrack*) tracks ->At(iTrack) ; | |
752 | trackParams = track ->GetTrackParamAtHit(); | |
753 | hits = track ->GetHitForRecAtHit() ; | |
754 | ||
755 | //Loop on hits | |
756 | nHits = (Int_t) hits -> GetEntriesFast(); | |
757 | for ( iHit = 0; iHit < nHits; ++iHit ) | |
758 | { | |
759 | trackParam = (AliMUONTrackParam*) trackParams -> At(iHit); | |
760 | hit = (AliMUONHitForRec*) hits -> At(iHit); | |
761 | chamberNbr = hit -> GetChamberNumber(); | |
762 | if ( chamberNbr >= firstChamber-1 ) { | |
763 | if ( chamberNbr < lastChamber ) { | |
764 | //Positions | |
765 | Double_t traX, traY; | |
766 | Double_t hitX, hitY; | |
767 | ||
768 | //Resolution | |
769 | Double_t deltaX, deltaY; | |
770 | ||
771 | traX = trackParam -> GetNonBendingCoor(); | |
772 | traY = trackParam -> GetBendingCoor(); | |
773 | hitX = hit -> GetNonBendingCoor(); | |
774 | hitY = hit -> GetBendingCoor(); | |
775 | ||
776 | deltaX = traX - hitX; | |
777 | deltaY = traY - hitY; | |
778 | ||
779 | hDeltax -> Fill (deltaX); | |
780 | hDeltay -> Fill (deltaY); | |
781 | hDelta3D -> Fill (deltaX, deltaY); | |
782 | } | |
783 | } | |
784 | } | |
785 | //End loop on hits | |
786 | } | |
787 | //End loop on tracks | |
788 | muondata.ResetRecTracks(); | |
789 | if ( event2Check!=0 ) | |
790 | iEvent=nEvents; | |
791 | } | |
792 | //End loop on events | |
793 | ||
794 | char hXtitle[80]; | |
795 | char hYtitle[80]; | |
796 | char h3title[80]; | |
797 | ||
798 | if ( firstChamber == lastChamber ) { | |
799 | sprintf(hXtitle,"DeltaX Chamber %d",firstChamber); | |
800 | sprintf(hYtitle,"DeltaY Chamber %d",lastChamber); | |
801 | sprintf(h3title,"DeltaX and Delta Y Chamber %d",lastChamber); | |
802 | } | |
803 | else{ | |
804 | sprintf(hXtitle,"DeltaX Chamber %d to %d",firstChamber,lastChamber); | |
805 | sprintf(hYtitle,"DeltaY Chamber %d to %d",firstChamber,lastChamber); | |
806 | sprintf(h3title,"DeltaX and Delta Y Chamber %d to %d",firstChamber,lastChamber); | |
807 | } | |
808 | ||
809 | ||
810 | hDeltax -> SetTitle (hXtitle); | |
811 | hDeltay -> SetTitle (hYtitle); | |
812 | hDelta3D -> SetTitle (h3title); | |
813 | ||
814 | Double_t rmsX = hDeltax -> GetRMS (); | |
815 | Double_t errX = hDeltax -> GetRMSError(); | |
816 | ||
817 | TF1 *fY = new TF1("fY","gaus",-0.3, 0.3); | |
818 | hDeltay -> Fit("fY","R","E"); | |
819 | Double_t sigY = fY -> GetParameter(2); | |
820 | Double_t errY = fY -> GetParError (2); | |
821 | ||
822 | if ( firstChamber == lastChamber ) { | |
823 | cout<<"\nThe resolution (Non Bending plane) of the chamber "<<firstChamber<<" = "<<rmsX<<" cm +-"<<errX; | |
824 | cout<<"\nThe resolution (Bending plane) of the chamber "<<firstChamber<<" = "<<sigY<<" cm +- "<<errY; | |
825 | } | |
826 | else { | |
827 | cout<<"\nThe resolution (Non Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<rmsX<<" cm +- "<<errX; | |
828 | cout<<"\nThe resolution (Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<sigY<<" cm +- "<<errY; | |
829 | } | |
830 | ||
831 | cout<<"\n\n\n"; | |
832 | MUONLoader->UnloadTracks(); | |
833 | ||
834 | } | |
835 | ||
836 | ||
837 | ||
838 | ||
839 | ||
840 | ||
841 | ||
842 | /*****************************************************************************************************************/ | |
843 | /*****************************************************************************************************************/ | |
844 | /*RESOLUTION IN FUNCTION OF PHI*/ | |
845 | ||
846 | void resolutionPhi( Int_t event2Check=0, char * filename="galice.root" ) | |
847 | { | |
848 | cout<<"\nChamber(s) chosen;\nFirst chamber:"; | |
849 | cin>>firstChamber; | |
850 | cout<<"Last chamber:"; | |
851 | cin>>lastChamber; | |
852 | cout<<"\n\n"; | |
853 | ||
854 | TH1F * hDeltaX[5]; | |
855 | TH1F * hDeltaY[5]; | |
856 | ||
857 | hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3); | |
858 | hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3); | |
859 | hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3); | |
860 | hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3); | |
861 | hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3); | |
862 | ||
863 | hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1); | |
864 | hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1); | |
865 | hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1); | |
866 | hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1); | |
867 | hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1); | |
868 | ||
869 | ||
870 | //Creating Run Loader and openning file containing Hits | |
871 | //-------------------------------------------------------------------------------------------------------------- | |
872 | AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ"); | |
873 | if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;} | |
874 | ||
875 | ||
876 | //Getting MUONLoader | |
877 | //-------------------------------------------------------------------------------------------------------------- | |
878 | AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader"); | |
879 | MUONLoader -> LoadTracks("READ"); | |
880 | MUONLoader -> LoadRecPoints("READ"); | |
881 | ||
882 | ||
883 | //Creating a MUON data container | |
884 | //-------------------------------------------------------------------------------------------------------------- | |
885 | AliMUONData muondata(MUONLoader,"MUON","MUON"); | |
886 | ||
887 | ||
888 | nEvents = (Int_t) RunLoader -> GetNumberOfEvents(); | |
889 | ||
890 | //Loop on events | |
891 | for ( iEvent = 0; iEvent < nEvents; ++iEvent ) | |
892 | { | |
893 | if ( event2Check!=0 ) | |
894 | iEvent=event2Check; | |
895 | printf("\r>>> Event %d ",iEvent); | |
896 | RunLoader->GetEvent(iEvent); | |
897 | ||
898 | //Addressing | |
899 | muondata.SetTreeAddress("RT"); | |
900 | muondata.GetRecTracks(); | |
901 | tracks = muondata.RecTracks(); | |
902 | ||
903 | //Loop on track | |
904 | nTracks = (Int_t) tracks -> GetEntriesFast(); | |
905 | for ( iTrack = 0; iTrack < nTracks; ++iTrack ) | |
906 | { | |
907 | track = (AliMUONTrack*) tracks ->At(iTrack) ; | |
908 | trackParams = track ->GetTrackParamAtHit(); | |
909 | hits = track ->GetHitForRecAtHit() ; | |
910 | ||
911 | //Loop on hits | |
912 | nHits = (Int_t) hits -> GetEntriesFast(); | |
913 | for ( iHit = 0; iHit < nHits; ++iHit ) | |
914 | { | |
915 | trackParam = (AliMUONTrackParam*) trackParams -> At(iHit); | |
916 | hit = (AliMUONHitForRec* ) hits -> At(iHit); | |
917 | chamberNbr = hit -> GetChamberNumber(); | |
918 | if ( chamberNbr >= firstChamber -1 ) { | |
919 | if ( chamberNbr < lastChamber ) { | |
920 | //Positions | |
921 | Double_t traX, traY; | |
922 | Double_t hitX, hitY; | |
923 | Double_t aveY, aveX; | |
924 | ||
925 | //Angles | |
926 | Double_t phi; | |
927 | Int_t iPhi; | |
928 | ||
929 | //Resolution | |
930 | Double_t deltaX; | |
931 | Double_t deltaY; | |
932 | ||
933 | traX = trackParam -> GetNonBendingCoor(); | |
934 | traY = trackParam -> GetBendingCoor() ; | |
935 | ||
936 | hitX = hit -> GetNonBendingCoor(); | |
937 | hitY = hit -> GetBendingCoor() ; | |
938 | ||
939 | aveX = 1./2.*(traX + hitX); | |
940 | aveY = 1./2.*(traY + hitY); | |
941 | ||
942 | //The calculation of phi: | |
943 | phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX )); | |
944 | ||
945 | if ( phi < 0.) phi = 360 - abs(phi); | |
946 | iPhi = int( phi/72. ); | |
947 | ||
948 | deltaX = traX - hitX; | |
949 | deltaY = traY - hitY; | |
950 | ||
951 | hDeltaX [iPhi] -> Fill( deltaX ); | |
952 | hDeltaY [iPhi] -> Fill( deltaY ); | |
953 | } | |
954 | } | |
955 | ||
956 | } | |
957 | //End loop on hits | |
958 | ||
959 | } | |
960 | //End loop on tracks | |
961 | muondata.ResetRecTracks(); | |
962 | if ( event2Check!=0 ) | |
963 | iEvent=nEvents; | |
964 | ||
965 | } | |
966 | //End loop on events | |
967 | ||
968 | ||
969 | Int_t iPhi; | |
970 | Int_t iPhiMax = 5; | |
971 | Int_t phiMin, phiMax; | |
972 | ||
973 | Float_t phi[5]; | |
974 | Float_t sigmaY[5]; | |
975 | Float_t errSY [5]; | |
976 | Float_t rmsX [5]; | |
977 | Float_t errSX [5]; | |
978 | Float_t errPh [5]; | |
979 | ||
980 | for ( iPhi = 0; iPhi < iPhiMax; ++iPhi ) | |
981 | { | |
982 | char hXtitle[80]; | |
983 | char hYtitle[80]; | |
984 | ||
985 | phiMin = 72*iPhi ; | |
986 | phiMax = 72*iPhi + 72; | |
987 | ||
988 | TF1 *fY = new TF1("fY","gaus",-0.3, 0.3); | |
989 | ||
990 | if ( firstChamber == lastChamber ) { | |
991 | sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d",phiMin,phiMax,firstChamber); | |
992 | sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d",phiMin,phiMax,lastChamber); | |
993 | } | |
994 | else{ | |
995 | sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber); | |
996 | sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber); | |
997 | } | |
998 | ||
999 | hDeltaY [iPhi] -> SetTitle(hYtitle); | |
1000 | hDeltaX [iPhi] -> SetTitle(hXtitle); | |
1001 | ||
1002 | hDeltaY [iPhi] -> Fit("fY","R","E"); | |
1003 | sigmaY [iPhi] = fY -> GetParameter(2) ; | |
1004 | errSY [iPhi] = fY -> GetParError(2) ; | |
1005 | ||
1006 | rmsX [iPhi] = hDeltaX [iPhi] -> GetRMS(); | |
1007 | errSX [iPhi] = hDeltaX [iPhi] -> GetRMSError(); | |
1008 | ||
1009 | phi [iPhi] = 72*iPhi + 36 ; | |
1010 | errPh [iPhi] = 36; | |
1011 | } | |
1012 | ||
1013 | //-------------------------------------------------------------------------------------------------------------- | |
1014 | //For plotting resolution in function of the angle of incidence | |
1015 | ||
1016 | TCanvas * c1 = new TCanvas("c1","",200,10,700,500); | |
1017 | c1-> Divide(1,2); | |
1018 | c1->cd(1); | |
1019 | ||
1020 | TGraphErrors * GraphX = new TGraphErrors( iPhiMax, phi, rmsX, errPh, errSX); | |
1021 | GraphX->SetTitle("Resolution en X (Phi)"); | |
1022 | GraphX->Draw("ALP"); | |
1023 | ||
1024 | c1->cd(2); | |
1025 | TGraphErrors * GraphY = new TGraphErrors( iPhiMax, phi, sigmaY, errPh, errSY); | |
1026 | GraphY->SetTitle("Resolution en Y (Phi)"); | |
1027 | GraphY->Draw("ALP"); | |
1028 | ||
1029 | cout<<"\n\n\n"; | |
1030 | ||
1031 | MUONLoader->UnloadTracks(); | |
1032 | ||
1033 | } | |
1034 | ||
1035 | ||
1036 | ||
1037 | ||
1038 | ||
1039 | ||
1040 | ||
1041 | /*****************************************************************************************************************/ | |
1042 | /*****************************************************************************************************************/ | |
1043 | /*RESOLUTION IN FUNCTION OF THETA*/ | |
1044 | ||
1045 | void resolutionTheta( Int_t event2Check=0, char * filename="galice.root" ) | |
1046 | { | |
1047 | cout<<"\nChamber(s) chosen;\nFirst chamber:"; | |
1048 | cin>>firstChamber; | |
1049 | cout<<"Last chamber:"; | |
1050 | cin>>lastChamber; | |
1051 | cout<<"\n\n"; | |
1052 | ||
1053 | TH1F * hDeltaX[9]; | |
1054 | TH1F * hDeltaY[9]; | |
1055 | ||
1056 | hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3); | |
1057 | hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3); | |
1058 | hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3); | |
1059 | hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3); | |
1060 | hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3); | |
1061 | hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3); | |
1062 | hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3); | |
1063 | hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3); | |
1064 | hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3); | |
1065 | ||
1066 | hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1); | |
1067 | hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1); | |
1068 | hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1); | |
1069 | hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1); | |
1070 | hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1); | |
1071 | hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1); | |
1072 | hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1); | |
1073 | hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1); | |
1074 | hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1); | |
1075 | ||
1076 | ||
1077 | //Creating Run Loader and openning file containing Hits | |
1078 | //-------------------------------------------------------------------------------------------------------------- | |
1079 | AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ"); | |
1080 | if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;} | |
1081 | ||
1082 | ||
1083 | //Getting MUONLoader | |
1084 | //-------------------------------------------------------------------------------------------------------------- | |
1085 | AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader"); | |
1086 | MUONLoader -> LoadTracks("READ"); | |
1087 | MUONLoader -> LoadRecPoints("READ"); | |
1088 | ||
1089 | ||
1090 | //Creating a MUON data container | |
1091 | //-------------------------------------------------------------------------------------------------------------- | |
1092 | AliMUONData muondata(MUONLoader,"MUON","MUON"); | |
1093 | ||
1094 | ||
1095 | nEvents = (Int_t) RunLoader -> GetNumberOfEvents(); | |
1096 | ||
1097 | //Loop on events | |
1098 | for ( iEvent = 0; iEvent < nEvents; ++iEvent ) | |
1099 | { | |
1100 | if ( event2Check!=0 ) | |
1101 | iEvent=event2Check; | |
1102 | printf("\r>>> Event %d ",iEvent); | |
1103 | RunLoader->GetEvent(iEvent); | |
1104 | ||
1105 | //Addressing | |
1106 | muondata.SetTreeAddress("RT"); | |
1107 | muondata.GetRecTracks(); | |
1108 | tracks = muondata.RecTracks(); | |
1109 | ||
1110 | //Loop on track | |
1111 | nTracks = (Int_t) tracks -> GetEntriesFast(); | |
1112 | for ( iTrack = 0; iTrack < nTracks; ++iTrack ) | |
1113 | { | |
1114 | track = (AliMUONTrack*) tracks ->At(iTrack) ; | |
1115 | trackParams = track ->GetTrackParamAtHit(); | |
1116 | hits = track ->GetHitForRecAtHit() ; | |
1117 | ||
1118 | //Loop on hits | |
1119 | nHits = (Int_t) hits -> GetEntriesFast(); | |
1120 | for ( iHit = 0; iHit < nHits; ++iHit ) | |
1121 | { | |
1122 | trackParam = (AliMUONTrackParam*) trackParams -> At(iHit); | |
1123 | hit = (AliMUONHitForRec* ) hits -> At(iHit); | |
1124 | chamberNbr = hit -> GetChamberNumber(); | |
1125 | if ( chamberNbr >= firstChamber -1 ) { | |
1126 | if ( chamberNbr < lastChamber ) { | |
1127 | //Positions | |
1128 | Double_t traX, traY; | |
1129 | Double_t hitX, hitY, hitZ; | |
1130 | Double_t aveY, aveX; | |
1131 | ||
1132 | //Angles | |
1133 | Double_t theta; | |
1134 | Int_t iTheta; | |
1135 | ||
1136 | //Resolution | |
1137 | Double_t deltaX; | |
1138 | Double_t deltaY; | |
1139 | ||
1140 | traX = trackParam -> GetNonBendingCoor(); | |
1141 | traY = trackParam -> GetBendingCoor() ; | |
1142 | ||
1143 | hitX = hit -> GetNonBendingCoor(); | |
1144 | hitY = hit -> GetBendingCoor() ; | |
1145 | hitZ = hit -> GetZ(); | |
1146 | ||
1147 | aveX = 1./2.*(traX + hitX); | |
1148 | aveY = 1./2.*(traY + hitY); | |
1149 | ||
1150 | //The calculation of theta, theta is in fact 180 - Theta(The polar angle): | |
1151 | theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ )); | |
1152 | ||
1153 | iTheta = int( theta ); | |
1154 | if( theta > 10 ) iTheta = 9; | |
1155 | if( theta < 1 ) iTheta = 1; | |
1156 | ||
1157 | deltaX = traX - hitX; | |
1158 | deltaY = traY - hitY; | |
1159 | ||
1160 | hDeltaX [iTheta-1] -> Fill( deltaX ); | |
1161 | hDeltaY [iTheta-1] -> Fill( deltaY ); | |
1162 | } | |
1163 | } | |
1164 | ||
1165 | } | |
1166 | //End loop on hits | |
1167 | ||
1168 | } | |
1169 | //End loop on tracks | |
1170 | muondata.ResetRecTracks(); | |
1171 | if ( event2Check!=0 ) | |
1172 | iEvent=nEvents; | |
1173 | ||
1174 | } | |
1175 | //End loop on events | |
1176 | ||
1177 | ||
1178 | Int_t iTheta; | |
1179 | Int_t iThetaMax = 9; | |
1180 | Int_t thetaMin, thetaMax; | |
1181 | ||
1182 | Float_t theta [9]; | |
1183 | Float_t sigmaY[9]; | |
1184 | Float_t errSY [9]; | |
1185 | Float_t rmsX [9]; | |
1186 | Float_t errSX [9]; | |
1187 | Float_t ErrTh [9]; | |
1188 | ||
1189 | for ( iTheta = 0; iTheta < iThetaMax; ++iTheta ) | |
1190 | { | |
1191 | char hXtitle[80]; | |
1192 | char hYtitle[80]; | |
1193 | ||
1194 | //To find the polar angle | |
1195 | thetaMin = 178 - iTheta; | |
1196 | thetaMax = 179 - iTheta; | |
1197 | ||
1198 | TF1 *fY = new TF1("fY","gaus",-0.3, 0.3); | |
1199 | ||
1200 | if ( firstChamber == lastChamber ) { | |
1201 | sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber); | |
1202 | sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber); | |
1203 | } | |
1204 | else{ | |
1205 | sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber); | |
1206 | sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber); | |
1207 | } | |
1208 | ||
1209 | hDeltaY [iTheta] -> SetTitle(hYtitle); | |
1210 | hDeltaX [iTheta] -> SetTitle(hXtitle); | |
1211 | ||
1212 | hDeltaY [iTheta] -> Fit("fY","R","E"); | |
1213 | sigmaY [iTheta] = fY -> GetParameter(2) ; | |
1214 | errSY [iTheta] = fY -> GetParError(2) ; | |
1215 | ||
1216 | rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS(); | |
1217 | errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError(); | |
1218 | ||
1219 | theta [iTheta] = 178.5 - iTheta; | |
1220 | ErrTh [iTheta] = 0.5; | |
1221 | } | |
1222 | ||
1223 | //-------------------------------------------------------------------------------------------------------------- | |
1224 | //For plotting resolution in function of the angle of incidence | |
1225 | ||
1226 | TCanvas * c1 = new TCanvas("c1","",200,10,700,500); | |
1227 | c1-> Divide(1,2); | |
1228 | c1->cd(1); | |
1229 | ||
1230 | TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, ErrTh, errSX); | |
1231 | GraphX->SetTitle("Resolution en X (Theta)"); | |
1232 | GraphX->Draw("ALP"); | |
1233 | ||
1234 | c1->cd(2); | |
1235 | TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, ErrTh, errSY); | |
1236 | GraphY->SetTitle("Resolution en Y (Theta)"); | |
1237 | GraphY->Draw("ALP"); | |
1238 | ||
1239 | cout<<"\n\n\n"; | |
1240 | MUONLoader->UnloadTracks(); | |
1241 | ||
1242 | } | |
1243 | ||
1244 | ||
1245 | ||
1246 | ||
1247 | ||
1248 | /*****************************************************************************************************************/ | |
1249 | /*****************************************************************************************************************/ | |
1250 | /*RESOLUTION IN FUNCTION OF THETA OF INCIDENCE*/ | |
1251 | ||
1252 | void resolutionThetaI( Int_t event2Check=0, char * filename="galice.root" ) | |
1253 | { | |
1254 | cout<<"\nChamber(s) chosen;\nFirst chamber:"; | |
1255 | cin>>firstChamber; | |
1256 | cout<<"Last chamber:"; | |
1257 | cin>>lastChamber; | |
1258 | cout<<"\n\n"; | |
1259 | ||
1260 | TH1F * hDeltaX[11]; | |
1261 | TH1F * hDeltaY[11]; | |
1262 | ||
1263 | hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3); | |
1264 | hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3); | |
1265 | hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3); | |
1266 | hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3); | |
1267 | hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3); | |
1268 | hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3); | |
1269 | hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3); | |
1270 | hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3); | |
1271 | hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3); | |
1272 | hDeltaX[9] = new TH1F ("hDeltaX9" , "No Interest", 100, -0.3, 0.3); | |
1273 | hDeltaX[10]= new TH1F ("hDeltaX10", "No Interest", 100, -0.3, 0.3); | |
1274 | ||
1275 | hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1); | |
1276 | hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1); | |
1277 | hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1); | |
1278 | hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1); | |
1279 | hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1); | |
1280 | hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1); | |
1281 | hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1); | |
1282 | hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1); | |
1283 | hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1); | |
1284 | hDeltaY[9] = new TH1F ("hDeltaY9" , "No Interest", 500, -0.1, 0.1); | |
1285 | hDeltaY[10]= new TH1F ("hDeltaY10", "No Interest", 500, -0.1, 0.1); | |
1286 | ||
1287 | ||
1288 | //Creating Run Loader and openning file containing Hits | |
1289 | //-------------------------------------------------------------------------------------------------------------- | |
1290 | AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ"); | |
1291 | if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;} | |
1292 | ||
1293 | ||
1294 | //Getting MUONLoader | |
1295 | //-------------------------------------------------------------------------------------------------------------- | |
1296 | AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader"); | |
1297 | MUONLoader -> LoadTracks("READ"); | |
1298 | MUONLoader -> LoadRecPoints("READ"); | |
1299 | ||
1300 | ||
1301 | //Creating a MUON data container | |
1302 | //-------------------------------------------------------------------------------------------------------------- | |
1303 | AliMUONData muondata(MUONLoader,"MUON","MUON"); | |
1304 | ||
1305 | ||
1306 | nEvents = (Int_t) RunLoader -> GetNumberOfEvents(); | |
1307 | ||
1308 | //Loop on events | |
1309 | for ( iEvent = 0; iEvent < nEvents; ++iEvent ) | |
1310 | { | |
1311 | if ( event2Check!=0 ) | |
1312 | iEvent=event2Check; | |
1313 | printf("\r>>> Event %d ",iEvent); | |
1314 | RunLoader->GetEvent(iEvent); | |
1315 | ||
1316 | //Addressing | |
1317 | muondata.SetTreeAddress("RT"); | |
1318 | muondata.GetRecTracks(); | |
1319 | tracks = muondata.RecTracks(); | |
1320 | ||
1321 | //Loop on track | |
1322 | nTracks = (Int_t) tracks -> GetEntriesFast(); | |
1323 | for ( iTrack = 0; iTrack < nTracks; ++iTrack ) | |
1324 | { | |
1325 | track = (AliMUONTrack*) tracks ->At(iTrack) ; | |
1326 | trackParams = track ->GetTrackParamAtHit(); | |
1327 | hits = track ->GetHitForRecAtHit() ; | |
1328 | ||
1329 | //Loop on hits | |
1330 | nHits = (Int_t) hits -> GetEntriesFast(); | |
1331 | for ( iHit = 0; iHit < nHits; ++iHit ) | |
1332 | { | |
1333 | trackParam = (AliMUONTrackParam*) trackParams -> At(iHit); | |
1334 | hit = (AliMUONHitForRec* ) hits -> At(iHit); | |
1335 | chamberNbr = hit -> GetChamberNumber(); | |
1336 | if ( chamberNbr >= firstChamber - 1 ) { | |
1337 | if ( chamberNbr < lastChamber ) { | |
1338 | //Momentum | |
1339 | Double_t Px, Py, Pz, Pr; | |
1340 | ||
1341 | //Positions | |
1342 | Double_t traX, traY; | |
1343 | Double_t hitX, hitY; | |
1344 | ||
1345 | //Resolution | |
1346 | Double_t deltaX; | |
1347 | Double_t deltaY; | |
1348 | ||
1349 | //Angle | |
1350 | Double_t theta; | |
1351 | Int_t iTheta; | |
1352 | ||
1353 | Px = trackParam -> Px() ; | |
1354 | Py = trackParam -> Py() ; | |
1355 | Pz = trackParam -> Pz() ; | |
1356 | Pr = sqrt( Px*Px + Py*Py ); | |
1357 | ||
1358 | traX = trackParam -> GetNonBendingCoor(); | |
1359 | traY = trackParam -> GetBendingCoor(); | |
1360 | ||
1361 | hitX = hit -> GetNonBendingCoor(); | |
1362 | hitY = hit -> GetBendingCoor(); | |
1363 | ||
1364 | //The calculation of theta, the angle of incidence: | |
1365 | theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr)) | |
1366 | ; | |
1367 | if ( theta < 79 ) iTheta = 0; | |
1368 | else iTheta = int( theta - 79.); | |
1369 | ||
1370 | deltaX = traX - hitX; | |
1371 | deltaY = traY - hitY; | |
1372 | ||
1373 | hDeltaX [iTheta] -> Fill (deltaX); | |
1374 | hDeltaY [iTheta] -> Fill (deltaY); | |
1375 | } | |
1376 | } | |
1377 | ||
1378 | } | |
1379 | //End loop on hits | |
1380 | ||
1381 | } | |
1382 | //End loop on tracks | |
1383 | muondata.ResetRecTracks(); | |
1384 | if ( event2Check!=0 ) | |
1385 | iEvent=nEvents; | |
1386 | ||
1387 | } | |
1388 | //End loop on events | |
1389 | ||
1390 | ||
1391 | Int_t iTheta; | |
1392 | Int_t iThetaMax = 11; | |
1393 | Int_t thetaMin, thetaMax; | |
1394 | ||
1395 | Float_t theta [11]; | |
1396 | Float_t sigmaY[11]; | |
1397 | Float_t errSY [11]; | |
1398 | Float_t rmsX [11]; | |
1399 | Float_t errSX [11]; | |
1400 | Float_t errTh [11]; | |
1401 | ||
1402 | for ( iTheta = 0; iTheta < iThetaMax; ++iTheta ) | |
1403 | { | |
1404 | char hXtitle[80]; | |
1405 | char hYtitle[80]; | |
1406 | ||
1407 | thetaMin = iTheta + 79; | |
1408 | thetaMax = iTheta + 80; | |
1409 | ||
1410 | TF1 *fY = new TF1("fY","gaus",-0.3, 0.3); | |
1411 | ||
1412 | if ( firstChamber == lastChamber ) { | |
1413 | sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber); | |
1414 | sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber); | |
1415 | } | |
1416 | else{ | |
1417 | sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber); | |
1418 | sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber); | |
1419 | } | |
1420 | ||
1421 | hDeltaY [iTheta] -> SetTitle(hYtitle); | |
1422 | hDeltaX [iTheta] -> SetTitle(hXtitle); | |
1423 | ||
1424 | hDeltaY [iTheta] -> Fit("fY","R","E"); | |
1425 | sigmaY [iTheta] = fY -> GetParameter(2) ; | |
1426 | errSY [iTheta] = fY -> GetParError(2) ; | |
1427 | ||
1428 | rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS(); | |
1429 | errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError(); | |
1430 | ||
1431 | theta [iTheta] = iTheta + 79.5 ; | |
1432 | errTh [iTheta] = 0.5; | |
1433 | } | |
1434 | ||
1435 | //-------------------------------------------------------------------------------------------------------------- | |
1436 | //For plotting resolution in function of the angle of incidence | |
1437 | ||
1438 | TCanvas * c1 = new TCanvas("c1","",200,10,700,500); | |
1439 | c1-> Divide(1,2); | |
1440 | c1->cd(1); | |
1441 | ||
1442 | TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, errTh, errSX); | |
1443 | GraphX->SetTitle("Resolution en X (Theta)"); | |
1444 | GraphX->Draw("ALP"); | |
1445 | ||
1446 | c1->cd(2); | |
1447 | TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, errTh, errSY); | |
1448 | GraphY->SetTitle("Resolution en Y (Theta)"); | |
1449 | GraphY->Draw("ALP"); | |
1450 | ||
1451 | cout<<"\n\n\n"; | |
1452 | MUONLoader->UnloadTracks(); | |
1453 | ||
1454 | } | |
1455 | ||
1456 | ||
1457 | ||
1458 |