]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/MUONResoEffChamber.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removing obsolete macro
[u/mrichter/AliRoot.git] / MUON / MUONResoEffChamber.C
CommitLineData
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
94void 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
200void 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
464void 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
697void 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
846void 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
1045void 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
1252void 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