1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
19 /// \file MUONResoEffChamber.C
20 /// \brief Macro to calculate the resolution and the efficiency of chamber(s)
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
28 /// \author Nicolas Le Bris, Subatech
30 #if !defined(__CINT__) || defined(__MAKECINT__)
34 #include "TClonesArray.h"
37 #include "TParticle.h"
48 #include "TGraphErrors.h"
56 #include "AliRunLoader.h"
57 #include "AliHeader.h"
58 #include "AliLoader.h"
59 #include "AliTracker.h"
66 #include "AliMUONData.h"
67 #include "AliMUONConstants.h"
69 #include "AliMUONHit.h"
70 #include "AliMUONHitForRec.h"
71 #include "AliMUONTrack.h"
72 #include "AliMUONTrackParam.h"
73 #include "AliMUONTrackExtrap.h"
75 #include "AliMpVSegmentation.h"
76 #include "AliMpIntPair.h"
77 #include "AliMpDEManager.h"
82 const Double_t kInvPi = 1./3.14159265;
88 Int_t nEvents, iEvent;
90 Int_t nTracks, iTrack;
94 Int_t firstChamber, lastChamber;
97 AliMUONTrack * track ;
98 AliMUONHitForRec * hit = 0;
99 AliMUONTrackParam * trackParam = 0;
101 TClonesArray * tracks ;
102 TClonesArray * trackParams;
103 TClonesArray * hits ;
110 /*****************************************************************************************************************/
111 /*****************************************************************************************************************/
114 void efficiency( Int_t event2Check=0, char * filename="galice.root" )
116 cout<<"\nChamber(s) chosen;\nFirst chamber:";
118 cout<<"Last chamber:";
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;}
129 //--------------------------------------------------------------------------------------------------------------
130 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
131 MUONLoader -> LoadTracks("READ");
132 MUONLoader -> LoadRecPoints("READ");
135 //Creating a MUON data container
136 //--------------------------------------------------------------------------------------------------------------
137 AliMUONData muondata(MUONLoader,"MUON","MUON");
140 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
144 Int_t chamber[10] = {0};
145 Int_t detEltNew, detElt;
148 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
150 if ( event2Check!=0 )
152 printf("\r>>> Event %d ",iEvent);
153 RunLoader -> GetEvent(iEvent);
155 muondata.SetTreeAddress("RT");
156 muondata.GetRecTracks();
157 tracks = muondata.RecTracks();
160 nTracks = (Int_t) tracks -> GetEntriesFast();
161 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
163 track = (AliMUONTrack*) tracks -> At(iTrack);
164 hits = track -> GetHitForRecAtHit();
167 nHits = (Int_t) hits -> GetEntriesFast();
169 for ( iHit = 0; iHit < nHits; ++iHit )
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 )
180 chamber[chamberNbr] = chamber[chamberNbr] + 1;
181 detEltOld = detEltNew;
190 muondata.ResetRecTracks();
191 if (event2Check != 0)
193 trackNb = trackNb + nTracks;
196 //--------------------------------------------------------------------------------------------------------------
199 for (Int_t i = firstChamber-1; i < lastChamber; ++i )
201 printf ("\nChamber %d has responded %d times on %d tracks", i+1, chamber[i], trackNb);
203 cout<<"\nEfficiency = ? (IS UNKNOWN)\n";
205 Double_t eff = chamber[i]*100./trackNb; cout<<"\nEfficiency = "<<eff<<" %\n";
209 MUONLoader->UnloadTracks();
216 /*****************************************************************************************************************/
217 /*****************************************************************************************************************/
218 /*EFFICIENCY IN FUNCTION OF THETA AND PHI*/
220 void efficiencyThetaPhi( Int_t event2Check=0, char * filename="galice.root" )
222 cout<<"\nChamber(s) chosen;\nFirst chamber:";
224 cout<<"Last chamber:";
228 Int_t eff [54] = {0};
229 Int_t trackN [54] = {0};
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;}
240 //--------------------------------------------------------------------------------------------------------------
241 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
242 MUONLoader -> LoadTracks("READ");
243 MUONLoader -> LoadRecPoints("READ");
246 //--------------------------------------------------------------------------------------------------------------
247 //Set mag field; waiting for mag field in CDB
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);
254 //--------------------------------------------------------------------------------------------------------------
255 //Set Field Map for track extrapolation
256 AliMUONTrackExtrap::SetField()
259 //Creating a MUON data container
260 //--------------------------------------------------------------------------------------------------------------
261 AliMUONData muondata(MUONLoader,"MUON","MUON");
264 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
267 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
269 if ( event2Check!=0 )
271 printf("\r>>> Event %d ",iEvent);
272 RunLoader->GetEvent(iEvent);
275 muondata.SetTreeAddress("RT");
276 muondata.GetRecTracks();
277 tracks = muondata.RecTracks();
280 nTracks = (Int_t) tracks -> GetEntriesFast();
281 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
283 track = (AliMUONTrack*) tracks ->At(iTrack) ;
284 trackParams = track ->GetTrackParamAtCluster();
285 hits = track ->GetHitForRecAtHit() ;
286 chamber = firstChamber-1;
291 nHits = (Int_t) hits->GetEntriesFast();
292 for ( iHit = 0; iHit<nHits; ++iHit )
294 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
295 hit = (AliMUONHitForRec* ) hits -> At(iHit);
296 chamberNbr = hit -> GetChamberNumber();
298 if ( chamberNbr == oldChamber )
301 oldChamber = chamberNbr;
302 if ( chamberNbr > chamber - k ) {
303 if ( chamber < lastChamber ) {
304 if ( chamberNbr == chamber ) {
306 Double_t traX, traY, traZ;
307 Double_t hitX, hitY, hitZ;
308 Double_t aveX, aveY ;
316 traX = trackParam -> GetNonBendingCoor();
317 traY = trackParam -> GetBendingCoor() ;
318 traZ = trackParam -> GetZ() ;
320 hitX = hit -> GetNonBendingCoor();
321 hitY = hit -> GetBendingCoor() ;
322 hitZ = hit -> GetZ() ;
324 aveX = 1./2.*(traX + hitX);
325 aveY = 1./2.*(traY + hitY);
327 //The calculation of phi:
328 phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));
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 ));
333 if ( phi < 0.) phi = 360 - fabs(phi);
334 iPhi = int( phi/72. );
335 iTheta = int( theta );
336 if( theta > 10 ) iTheta = 9;
337 if( theta < 1 ) iTheta = 1;
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;
346 Double_t chamberZpos;
347 Double_t exXpos, exYpos;
355 chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
356 AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
357 exXpos = (Double_t) trackParam->GetNonBendingCoor();
358 exYpos = (Double_t) trackParam->GetBendingCoor();
360 //The calculation of phi:
361 phi = 180. * kInvPi * (TMath::ATan2( exYpos, exXpos ));
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 ));
366 if ( phi < 0.) phi = 360. - fabs(phi);
367 iPhi = int( phi/72. );
368 iTheta = int( theta );
369 if( theta > 10 ) iTheta = 9;
370 if( theta < 1 ) iTheta = 1;
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;
381 if ( iHit == nHits-1 ) {
382 if ( chamber == lastChamber )
396 muondata.ResetRecTracks();
397 if ( event2Check!=0 )
403 TCanvas * c1 = new TCanvas();
404 TGraph2D* effPhiTheta = new TGraph2D();
405 Double_t efficiency = 0;
407 if ( firstChamber == lastChamber )
409 for ( Int_t ph = 0; ph < 5; ++ph )
411 for ( Int_t th = 1; th < 10; ++th )
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";
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;
422 if ( trackN[i] == 0 ) {
424 cout<<"Efficiency = ? % ( IS UNKNOWN )\n";
427 efficiency = 100.*e/n;
428 cout<<"Efficiency = "<<efficiency<<" %\n";
431 effPhiTheta -> SetPoint( i, p, t, efficiency);
437 for ( Int_t ph = 0; ph < 5; ++ph )
439 for ( Int_t th = 1; th < 10; ++th )
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";
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;
450 if ( trackN[i] == 0 ) {
452 cout<<"Efficiency = ? % ( IS UNKNOWN )\n";
455 efficiency = 100.*e/n;
456 cout<<"Efficiency = "<<efficiency<<" %\n";
459 effPhiTheta -> SetPoint( i, p, t, efficiency);
464 gStyle->SetPalette(1);
465 effPhiTheta -> Draw("surf1");
468 MUONLoader->UnloadTracks();
479 /*****************************************************************************************************************/
480 /*****************************************************************************************************************/
481 /*EFFICIENCY IN FUNCTION OF THETA OF INCIDENCE*/
484 void efficiencyThetaI( Int_t event2Check=0, char * filename="galice.root" )
486 cout<<"\nChamber(s) chosen;\nFirst chamber:";
488 cout<<"Last chamber:";
492 Int_t eff [12] = {0};
493 Int_t trackN [12] = {0};
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;}
505 //--------------------------------------------------------------------------------------------------------------
506 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
507 MUONLoader -> LoadTracks("READ");
508 MUONLoader -> LoadRecPoints("READ");
511 //--------------------------------------------------------------------------------------------------------------
512 //Set mag field; waiting for mag field in CDB
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");
520 //--------------------------------------------------------------------------------------------------------------
521 //Set Field Map for track extrapolation
522 AliMUONTrackExtrap::SetField();
525 //Creating a MUON data container
526 //--------------------------------------------------------------------------------------------------------------
527 AliMUONData muondata(MUONLoader,"MUON","MUON");
530 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
533 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
535 if ( event2Check!=0 )
537 printf("\r>>> Event %d ",iEvent);
538 RunLoader->GetEvent(iEvent);
541 muondata.SetTreeAddress("RT");
542 muondata.GetRecTracks();
543 tracks = muondata.RecTracks();
546 nTracks = (Int_t) tracks -> GetEntriesFast();
547 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
549 track = (AliMUONTrack*) tracks ->At(iTrack) ;
550 trackParams = track ->GetTrackParamAtCluster();
551 hits = track ->GetHitForRecAtHit() ;
552 chamber = firstChamber - 1;
557 nHits = (Int_t) hits -> GetEntriesFast();
558 for ( iHit = 0; iHit < nHits; ++iHit )
560 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
561 hit = (AliMUONHitForRec*) hits -> At(iHit);
562 chamberNbr = hit -> GetChamberNumber();
564 if ( chamberNbr == oldChamber )
567 oldChamber = chamberNbr;
568 if ( chamberNbr > chamber - k ) {
569 if ( chamber < lastChamber ) {
570 if ( chamberNbr == chamber ) {
572 Double_t Px, Py, Pz, Pr;
578 Px = trackParam -> Px() ;
579 Py = trackParam -> Py() ;
580 Pz = trackParam -> Pz() ;
581 Pr = sqrt( Px*Px + Py*Py );
583 //The calculation of theta, the angle of incidence:
584 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
586 if ( theta < 79 ) iTheta = 11;
587 else if ( theta < 90 ) iTheta = int( theta - 79.);
590 eff [iTheta] = eff [iTheta] + 1;
591 trackN [iTheta] = trackN [iTheta] + 1;
592 chamber = chamber + 1;
597 Double_t chamberZpos;
600 Double_t Px, Py, Pz, Pr;
606 chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
607 AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
609 Px = trackParam -> Px() ;
610 Py = trackParam -> Py() ;
611 Pz = trackParam -> Pz() ;
612 Pr = sqrt( Px*Px + Py*Py );
614 //The calculation of thetaI, the angle of incidence:
615 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
617 if ( theta < 79 ) iTheta = 11;
618 else if ( theta < 90 ) iTheta = int( theta - 79.);
621 eff [iTheta] = eff [iTheta] + 0;
622 trackN [iTheta] = trackN [iTheta] + 1;
623 chamber = chamber + 1;
630 if ( iHit == nHits-1 ) {
631 if ( chamber == lastChamber )
645 muondata.ResetRecTracks();
646 if ( event2Check!=0 )
652 Double_t efficiency [11];
656 TGraph * effTheta = new TGraph ();
658 if ( firstChamber == lastChamber ) {
659 for ( th = 0; th < 11; ++th )
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";
665 Double_t e = eff [th];
666 Double_t n = trackN [th];
669 efficiency [th] = 0.;
670 cout<<"Efficiency = ? % (IS UNKNOWN) \n";
673 efficiency [th] = 100.*e/n;
674 cout<<"Efficiency = "<<efficiency [th]<<" %\n";
680 for ( th = 0; th < 11; ++th )
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";
686 Double_t e = eff [th];
687 Double_t n = trackN [th];
690 efficiency [th] = 0.;
691 cout<<"Efficiency = ? % (IS UNKNOWN) \n";
694 efficiency [th] = 100.*e/n;
695 cout<<"Efficiency = "<<efficiency [th]<<" %\n";
700 TCanvas * c1 = new TCanvas ();
701 effTheta = new TGraph( i, t, efficiency );
703 effTheta -> Draw("ALP");
705 MUONLoader->UnloadTracks();
714 /*****************************************************************************************************************/
715 /*****************************************************************************************************************/
718 void resolution( Int_t event2Check=0, char * filename="galice.root" )
720 cout<<"\nChamber(s) chosen;\nFirst chamber:";
722 cout<<"Last chamber:";
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);
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;}
742 //--------------------------------------------------------------------------------------------------------------
743 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
744 MUONLoader -> LoadTracks("READ");
745 MUONLoader -> LoadRecPoints("READ");
748 //Creating a MUON data container
749 //--------------------------------------------------------------------------------------------------------------
750 AliMUONData muondata(MUONLoader,"MUON","MUON");
753 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
756 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
760 printf("\r>>> Event %d ",iEvent);
761 RunLoader->GetEvent(iEvent);
764 muondata.SetTreeAddress("RT");
765 muondata.GetRecTracks();
766 tracks = muondata.RecTracks();
769 nTracks = (Int_t) tracks -> GetEntriesFast();
770 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
772 track = (AliMUONTrack*) tracks ->At(iTrack) ;
773 trackParams = track ->GetTrackParamAtCluster();
774 hits = track ->GetHitForRecAtHit() ;
777 nHits = (Int_t) hits -> GetEntriesFast();
778 for ( iHit = 0; iHit < nHits; ++iHit )
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 ) {
790 Double_t deltaX, deltaY;
792 traX = trackParam -> GetNonBendingCoor();
793 traY = trackParam -> GetBendingCoor();
794 hitX = hit -> GetNonBendingCoor();
795 hitY = hit -> GetBendingCoor();
797 deltaX = traX - hitX;
798 deltaY = traY - hitY;
800 hDeltax -> Fill (deltaX);
801 hDeltay -> Fill (deltaY);
802 hDelta3D -> Fill (deltaX, deltaY);
809 muondata.ResetRecTracks();
810 if ( event2Check!=0 )
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);
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);
831 hDeltax -> SetTitle (hXtitle);
832 hDeltay -> SetTitle (hYtitle);
833 hDelta3D -> SetTitle (h3title);
835 Double_t rmsX = hDeltax -> GetRMS ();
836 Double_t errX = hDeltax -> GetRMSError();
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);
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;
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;
853 MUONLoader->UnloadTracks();
863 /*****************************************************************************************************************/
864 /*****************************************************************************************************************/
865 /*RESOLUTION IN FUNCTION OF PHI*/
867 void resolutionPhi( Int_t event2Check=0, char * filename="galice.root" )
869 cout<<"\nChamber(s) chosen;\nFirst chamber:";
871 cout<<"Last chamber:";
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);
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);
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;}
898 //--------------------------------------------------------------------------------------------------------------
899 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
900 MUONLoader -> LoadTracks("READ");
901 MUONLoader -> LoadRecPoints("READ");
904 //Creating a MUON data container
905 //--------------------------------------------------------------------------------------------------------------
906 AliMUONData muondata(MUONLoader,"MUON","MUON");
909 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
912 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
914 if ( event2Check!=0 )
916 printf("\r>>> Event %d ",iEvent);
917 RunLoader->GetEvent(iEvent);
920 muondata.SetTreeAddress("RT");
921 muondata.GetRecTracks();
922 tracks = muondata.RecTracks();
925 nTracks = (Int_t) tracks -> GetEntriesFast();
926 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
928 track = (AliMUONTrack*) tracks ->At(iTrack) ;
929 trackParams = track ->GetTrackParamAtCluster();
930 hits = track ->GetHitForRecAtHit() ;
933 nHits = (Int_t) hits -> GetEntriesFast();
934 for ( iHit = 0; iHit < nHits; ++iHit )
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 ) {
954 traX = trackParam -> GetNonBendingCoor();
955 traY = trackParam -> GetBendingCoor() ;
957 hitX = hit -> GetNonBendingCoor();
958 hitY = hit -> GetBendingCoor() ;
960 aveX = 1./2.*(traX + hitX);
961 aveY = 1./2.*(traY + hitY);
963 //The calculation of phi:
964 phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));
966 if ( phi < 0.) phi = 360. - fabs(phi);
967 iPhi = int( phi/72. );
969 deltaX = traX - hitX;
970 deltaY = traY - hitY;
972 hDeltaX [iPhi] -> Fill( deltaX );
973 hDeltaY [iPhi] -> Fill( deltaY );
982 muondata.ResetRecTracks();
983 if ( event2Check!=0 )
992 Int_t phiMin, phiMax;
1001 for ( iPhi = 0; iPhi < iPhiMax; ++iPhi )
1007 phiMax = 72*iPhi + 72;
1009 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
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);
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);
1020 hDeltaY [iPhi] -> SetTitle(hYtitle);
1021 hDeltaX [iPhi] -> SetTitle(hXtitle);
1023 hDeltaY [iPhi] -> Fit("fY","R","E");
1024 sigmaY [iPhi] = fY -> GetParameter(2) ;
1025 errSY [iPhi] = fY -> GetParError(2) ;
1027 rmsX [iPhi] = hDeltaX [iPhi] -> GetRMS();
1028 errSX [iPhi] = hDeltaX [iPhi] -> GetRMSError();
1030 phi [iPhi] = 72*iPhi + 36 ;
1034 //--------------------------------------------------------------------------------------------------------------
1035 //For plotting resolution in function of the angle of incidence
1037 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1041 TGraphErrors * GraphX = new TGraphErrors( iPhiMax, phi, rmsX, errPh, errSX);
1042 GraphX->SetTitle("Resolution en X (Phi)");
1043 GraphX->Draw("ALP");
1046 TGraphErrors * GraphY = new TGraphErrors( iPhiMax, phi, sigmaY, errPh, errSY);
1047 GraphY->SetTitle("Resolution en Y (Phi)");
1048 GraphY->Draw("ALP");
1052 MUONLoader->UnloadTracks();
1062 /*****************************************************************************************************************/
1063 /*****************************************************************************************************************/
1064 /*RESOLUTION IN FUNCTION OF THETA*/
1066 void resolutionTheta( Int_t event2Check=0, char * filename="galice.root" )
1068 cout<<"\nChamber(s) chosen;\nFirst chamber:";
1070 cout<<"Last chamber:";
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);
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);
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;}
1104 //Getting MUONLoader
1105 //--------------------------------------------------------------------------------------------------------------
1106 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1107 MUONLoader -> LoadTracks("READ");
1108 MUONLoader -> LoadRecPoints("READ");
1111 //Creating a MUON data container
1112 //--------------------------------------------------------------------------------------------------------------
1113 AliMUONData muondata(MUONLoader,"MUON","MUON");
1116 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1119 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1121 if ( event2Check!=0 )
1123 printf("\r>>> Event %d ",iEvent);
1124 RunLoader->GetEvent(iEvent);
1127 muondata.SetTreeAddress("RT");
1128 muondata.GetRecTracks();
1129 tracks = muondata.RecTracks();
1132 nTracks = (Int_t) tracks -> GetEntriesFast();
1133 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1135 track = (AliMUONTrack*) tracks ->At(iTrack) ;
1136 trackParams = track ->GetTrackParamAtCluster();
1137 hits = track ->GetHitForRecAtHit() ;
1140 nHits = (Int_t) hits -> GetEntriesFast();
1141 for ( iHit = 0; iHit < nHits; ++iHit )
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 ) {
1149 Double_t traX, traY;
1150 Double_t hitX, hitY, hitZ;
1151 Double_t aveY, aveX;
1161 traX = trackParam -> GetNonBendingCoor();
1162 traY = trackParam -> GetBendingCoor() ;
1164 hitX = hit -> GetNonBendingCoor();
1165 hitY = hit -> GetBendingCoor() ;
1166 hitZ = hit -> GetZ();
1168 aveX = 1./2.*(traX + hitX);
1169 aveY = 1./2.*(traY + hitY);
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 ));
1174 iTheta = int( theta );
1175 if( theta > 10 ) iTheta = 9;
1176 if( theta < 1 ) iTheta = 1;
1178 deltaX = traX - hitX;
1179 deltaY = traY - hitY;
1181 hDeltaX [iTheta-1] -> Fill( deltaX );
1182 hDeltaY [iTheta-1] -> Fill( deltaY );
1190 //End loop on tracks
1191 muondata.ResetRecTracks();
1192 if ( event2Check!=0 )
1196 //End loop on events
1200 Int_t iThetaMax = 9;
1201 Int_t thetaMin, thetaMax;
1210 for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1215 //To find the polar angle
1216 thetaMin = 178 - iTheta;
1217 thetaMax = 179 - iTheta;
1219 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
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);
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);
1230 hDeltaY [iTheta] -> SetTitle(hYtitle);
1231 hDeltaX [iTheta] -> SetTitle(hXtitle);
1233 hDeltaY [iTheta] -> Fit("fY","R","E");
1234 sigmaY [iTheta] = fY -> GetParameter(2) ;
1235 errSY [iTheta] = fY -> GetParError(2) ;
1237 rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS();
1238 errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1240 theta [iTheta] = 178.5 - iTheta;
1241 ErrTh [iTheta] = 0.5;
1244 //--------------------------------------------------------------------------------------------------------------
1245 //For plotting resolution in function of the angle of incidence
1247 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1251 TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, ErrTh, errSX);
1252 GraphX->SetTitle("Resolution en X (Theta)");
1253 GraphX->Draw("ALP");
1256 TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, ErrTh, errSY);
1257 GraphY->SetTitle("Resolution en Y (Theta)");
1258 GraphY->Draw("ALP");
1261 MUONLoader->UnloadTracks();
1269 /*****************************************************************************************************************/
1270 /*****************************************************************************************************************/
1271 /*RESOLUTION IN FUNCTION OF THETA OF INCIDENCE*/
1273 void resolutionThetaI( Int_t event2Check=0, char * filename="galice.root" )
1275 cout<<"\nChamber(s) chosen;\nFirst chamber:";
1277 cout<<"Last chamber:";
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);
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);
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;}
1315 //Getting MUONLoader
1316 //--------------------------------------------------------------------------------------------------------------
1317 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1318 MUONLoader -> LoadTracks("READ");
1319 MUONLoader -> LoadRecPoints("READ");
1322 //Creating a MUON data container
1323 //--------------------------------------------------------------------------------------------------------------
1324 AliMUONData muondata(MUONLoader,"MUON","MUON");
1327 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1330 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1332 if ( event2Check!=0 )
1334 printf("\r>>> Event %d ",iEvent);
1335 RunLoader->GetEvent(iEvent);
1338 muondata.SetTreeAddress("RT");
1339 muondata.GetRecTracks();
1340 tracks = muondata.RecTracks();
1343 nTracks = (Int_t) tracks -> GetEntriesFast();
1344 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1346 track = (AliMUONTrack*) tracks ->At(iTrack) ;
1347 trackParams = track ->GetTrackParamAtCluster();
1348 hits = track ->GetHitForRecAtHit() ;
1351 nHits = (Int_t) hits -> GetEntriesFast();
1352 for ( iHit = 0; iHit < nHits; ++iHit )
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 ) {
1360 Double_t Px, Py, Pz, Pr;
1363 Double_t traX, traY;
1364 Double_t hitX, hitY;
1374 Px = trackParam -> Px() ;
1375 Py = trackParam -> Py() ;
1376 Pz = trackParam -> Pz() ;
1377 Pr = sqrt( Px*Px + Py*Py );
1379 traX = trackParam -> GetNonBendingCoor();
1380 traY = trackParam -> GetBendingCoor();
1382 hitX = hit -> GetNonBendingCoor();
1383 hitY = hit -> GetBendingCoor();
1385 //The calculation of theta, the angle of incidence:
1386 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr))
1388 if ( theta < 79 ) iTheta = 0;
1389 else iTheta = int( theta - 79.);
1391 deltaX = traX - hitX;
1392 deltaY = traY - hitY;
1394 hDeltaX [iTheta] -> Fill (deltaX);
1395 hDeltaY [iTheta] -> Fill (deltaY);
1403 //End loop on tracks
1404 muondata.ResetRecTracks();
1405 if ( event2Check!=0 )
1409 //End loop on events
1413 Int_t iThetaMax = 11;
1414 Int_t thetaMin, thetaMax;
1423 for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1428 thetaMin = iTheta + 79;
1429 thetaMax = iTheta + 80;
1431 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
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);
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);
1442 hDeltaY [iTheta] -> SetTitle(hYtitle);
1443 hDeltaX [iTheta] -> SetTitle(hXtitle);
1445 hDeltaY [iTheta] -> Fit("fY","R","E");
1446 sigmaY [iTheta] = fY -> GetParameter(2) ;
1447 errSY [iTheta] = fY -> GetParError(2) ;
1449 rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS();
1450 errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1452 theta [iTheta] = iTheta + 79.5 ;
1453 errTh [iTheta] = 0.5;
1456 //--------------------------------------------------------------------------------------------------------------
1457 //For plotting resolution in function of the angle of incidence
1459 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1463 TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, errTh, errSX);
1464 GraphX->SetTitle("Resolution en X (Theta)");
1465 GraphX->Draw("ALP");
1468 TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, errTh, errSY);
1469 GraphY->SetTitle("Resolution en Y (Theta)");
1470 GraphY->Draw("ALP");
1473 MUONLoader->UnloadTracks();