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