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