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"
61 #include "AliMagFMaps.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 printf("Loading field map...\n");
249 if (!AliTracker::GetFieldMap()) {
250 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
251 AliTracker::SetFieldMap(field, kFALSE);}
254 //--------------------------------------------------------------------------------------------------------------
255 //Set Field Map for track extrapolation
256 AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
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 printf("Loading field map...\n");
514 if (!AliTracker::GetFieldMap()) {
515 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
516 AliTracker::SetFieldMap(field, kFALSE);}
519 //--------------------------------------------------------------------------------------------------------------
520 //Set Field Map for track extrapolation
521 AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
524 //Creating a MUON data container
525 //--------------------------------------------------------------------------------------------------------------
526 AliMUONData muondata(MUONLoader,"MUON","MUON");
529 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
532 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
534 if ( event2Check!=0 )
536 printf("\r>>> Event %d ",iEvent);
537 RunLoader->GetEvent(iEvent);
540 muondata.SetTreeAddress("RT");
541 muondata.GetRecTracks();
542 tracks = muondata.RecTracks();
545 nTracks = (Int_t) tracks -> GetEntriesFast();
546 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
548 track = (AliMUONTrack*) tracks ->At(iTrack) ;
549 trackParams = track ->GetTrackParamAtCluster();
550 hits = track ->GetHitForRecAtHit() ;
551 chamber = firstChamber - 1;
556 nHits = (Int_t) hits -> GetEntriesFast();
557 for ( iHit = 0; iHit < nHits; ++iHit )
559 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
560 hit = (AliMUONHitForRec*) hits -> At(iHit);
561 chamberNbr = hit -> GetChamberNumber();
563 if ( chamberNbr == oldChamber )
566 oldChamber = chamberNbr;
567 if ( chamberNbr > chamber - k ) {
568 if ( chamber < lastChamber ) {
569 if ( chamberNbr == chamber ) {
571 Double_t Px, Py, Pz, Pr;
577 Px = trackParam -> Px() ;
578 Py = trackParam -> Py() ;
579 Pz = trackParam -> Pz() ;
580 Pr = sqrt( Px*Px + Py*Py );
582 //The calculation of theta, the angle of incidence:
583 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
585 if ( theta < 79 ) iTheta = 11;
586 else if ( theta < 90 ) iTheta = int( theta - 79.);
589 eff [iTheta] = eff [iTheta] + 1;
590 trackN [iTheta] = trackN [iTheta] + 1;
591 chamber = chamber + 1;
596 Double_t chamberZpos;
599 Double_t Px, Py, Pz, Pr;
605 chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
606 AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
608 Px = trackParam -> Px() ;
609 Py = trackParam -> Py() ;
610 Pz = trackParam -> Pz() ;
611 Pr = sqrt( Px*Px + Py*Py );
613 //The calculation of thetaI, the angle of incidence:
614 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
616 if ( theta < 79 ) iTheta = 11;
617 else if ( theta < 90 ) iTheta = int( theta - 79.);
620 eff [iTheta] = eff [iTheta] + 0;
621 trackN [iTheta] = trackN [iTheta] + 1;
622 chamber = chamber + 1;
629 if ( iHit == nHits-1 ) {
630 if ( chamber == lastChamber )
644 muondata.ResetRecTracks();
645 if ( event2Check!=0 )
651 Double_t efficiency [11];
655 TGraph * effTheta = new TGraph ();
657 if ( firstChamber == lastChamber ) {
658 for ( th = 0; th < 11; ++th )
660 cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
661 cout<<"The chamber "<<firstChamber<<" has responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
664 Double_t e = eff [th];
665 Double_t n = trackN [th];
668 efficiency [th] = 0.;
669 cout<<"Efficiency = ? % (IS UNKNOWN) \n";
672 efficiency [th] = 100.*e/n;
673 cout<<"Efficiency = "<<efficiency [th]<<" %\n";
679 for ( th = 0; th < 11; ++th )
681 cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
682 cout<<"The chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
685 Double_t e = eff [th];
686 Double_t n = trackN [th];
689 efficiency [th] = 0.;
690 cout<<"Efficiency = ? % (IS UNKNOWN) \n";
693 efficiency [th] = 100.*e/n;
694 cout<<"Efficiency = "<<efficiency [th]<<" %\n";
699 TCanvas * c1 = new TCanvas ();
700 effTheta = new TGraph( i, t, efficiency );
702 effTheta -> Draw("ALP");
704 MUONLoader->UnloadTracks();
713 /*****************************************************************************************************************/
714 /*****************************************************************************************************************/
717 void resolution( Int_t event2Check=0, char * filename="galice.root" )
719 cout<<"\nChamber(s) chosen;\nFirst chamber:";
721 cout<<"Last chamber:";
729 hDeltax = new TH1F ("hDeltax " , "No Interest", 100, -0.3, 0.3);
730 hDeltay = new TH1F ("hDeltay " , "No Interest", 500, -0.1, 0.1);
731 hDelta3D = new TH2F ("hDelta3D" , "No Interest", 100, -0.3, 0.3, 500, -0.1, 0.1);
734 //Creating Run Loader and openning file containing Hits
735 //--------------------------------------------------------------------------------------------------------------
736 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
737 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
741 //--------------------------------------------------------------------------------------------------------------
742 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
743 MUONLoader -> LoadTracks("READ");
744 MUONLoader -> LoadRecPoints("READ");
747 //Creating a MUON data container
748 //--------------------------------------------------------------------------------------------------------------
749 AliMUONData muondata(MUONLoader,"MUON","MUON");
752 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
755 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
759 printf("\r>>> Event %d ",iEvent);
760 RunLoader->GetEvent(iEvent);
763 muondata.SetTreeAddress("RT");
764 muondata.GetRecTracks();
765 tracks = muondata.RecTracks();
768 nTracks = (Int_t) tracks -> GetEntriesFast();
769 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
771 track = (AliMUONTrack*) tracks ->At(iTrack) ;
772 trackParams = track ->GetTrackParamAtCluster();
773 hits = track ->GetHitForRecAtHit() ;
776 nHits = (Int_t) hits -> GetEntriesFast();
777 for ( iHit = 0; iHit < nHits; ++iHit )
779 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
780 hit = (AliMUONHitForRec*) hits -> At(iHit);
781 chamberNbr = hit -> GetChamberNumber();
782 if ( chamberNbr >= firstChamber-1 ) {
783 if ( chamberNbr < lastChamber ) {
789 Double_t deltaX, deltaY;
791 traX = trackParam -> GetNonBendingCoor();
792 traY = trackParam -> GetBendingCoor();
793 hitX = hit -> GetNonBendingCoor();
794 hitY = hit -> GetBendingCoor();
796 deltaX = traX - hitX;
797 deltaY = traY - hitY;
799 hDeltax -> Fill (deltaX);
800 hDeltay -> Fill (deltaY);
801 hDelta3D -> Fill (deltaX, deltaY);
808 muondata.ResetRecTracks();
809 if ( event2Check!=0 )
818 if ( firstChamber == lastChamber ) {
819 sprintf(hXtitle,"DeltaX Chamber %d",firstChamber);
820 sprintf(hYtitle,"DeltaY Chamber %d",lastChamber);
821 sprintf(h3title,"DeltaX and Delta Y Chamber %d",lastChamber);
824 sprintf(hXtitle,"DeltaX Chamber %d to %d",firstChamber,lastChamber);
825 sprintf(hYtitle,"DeltaY Chamber %d to %d",firstChamber,lastChamber);
826 sprintf(h3title,"DeltaX and Delta Y Chamber %d to %d",firstChamber,lastChamber);
830 hDeltax -> SetTitle (hXtitle);
831 hDeltay -> SetTitle (hYtitle);
832 hDelta3D -> SetTitle (h3title);
834 Double_t rmsX = hDeltax -> GetRMS ();
835 Double_t errX = hDeltax -> GetRMSError();
837 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
838 hDeltay -> Fit("fY","R","E");
839 Double_t sigY = fY -> GetParameter(2);
840 Double_t errY = fY -> GetParError (2);
842 if ( firstChamber == lastChamber ) {
843 cout<<"\nThe resolution (Non Bending plane) of the chamber "<<firstChamber<<" = "<<rmsX<<" cm +-"<<errX;
844 cout<<"\nThe resolution (Bending plane) of the chamber "<<firstChamber<<" = "<<sigY<<" cm +- "<<errY;
847 cout<<"\nThe resolution (Non Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<rmsX<<" cm +- "<<errX;
848 cout<<"\nThe resolution (Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<sigY<<" cm +- "<<errY;
852 MUONLoader->UnloadTracks();
862 /*****************************************************************************************************************/
863 /*****************************************************************************************************************/
864 /*RESOLUTION IN FUNCTION OF PHI*/
866 void resolutionPhi( Int_t event2Check=0, char * filename="galice.root" )
868 cout<<"\nChamber(s) chosen;\nFirst chamber:";
870 cout<<"Last chamber:";
877 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
878 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
879 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
880 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
881 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
883 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
884 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
885 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
886 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
887 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
890 //Creating Run Loader and openning file containing Hits
891 //--------------------------------------------------------------------------------------------------------------
892 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
893 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
897 //--------------------------------------------------------------------------------------------------------------
898 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
899 MUONLoader -> LoadTracks("READ");
900 MUONLoader -> LoadRecPoints("READ");
903 //Creating a MUON data container
904 //--------------------------------------------------------------------------------------------------------------
905 AliMUONData muondata(MUONLoader,"MUON","MUON");
908 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
911 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
913 if ( event2Check!=0 )
915 printf("\r>>> Event %d ",iEvent);
916 RunLoader->GetEvent(iEvent);
919 muondata.SetTreeAddress("RT");
920 muondata.GetRecTracks();
921 tracks = muondata.RecTracks();
924 nTracks = (Int_t) tracks -> GetEntriesFast();
925 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
927 track = (AliMUONTrack*) tracks ->At(iTrack) ;
928 trackParams = track ->GetTrackParamAtCluster();
929 hits = track ->GetHitForRecAtHit() ;
932 nHits = (Int_t) hits -> GetEntriesFast();
933 for ( iHit = 0; iHit < nHits; ++iHit )
935 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
936 hit = (AliMUONHitForRec* ) hits -> At(iHit);
937 chamberNbr = hit -> GetChamberNumber();
938 if ( chamberNbr >= firstChamber -1 ) {
939 if ( chamberNbr < lastChamber ) {
953 traX = trackParam -> GetNonBendingCoor();
954 traY = trackParam -> GetBendingCoor() ;
956 hitX = hit -> GetNonBendingCoor();
957 hitY = hit -> GetBendingCoor() ;
959 aveX = 1./2.*(traX + hitX);
960 aveY = 1./2.*(traY + hitY);
962 //The calculation of phi:
963 phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));
965 if ( phi < 0.) phi = 360. - fabs(phi);
966 iPhi = int( phi/72. );
968 deltaX = traX - hitX;
969 deltaY = traY - hitY;
971 hDeltaX [iPhi] -> Fill( deltaX );
972 hDeltaY [iPhi] -> Fill( deltaY );
981 muondata.ResetRecTracks();
982 if ( event2Check!=0 )
991 Int_t phiMin, phiMax;
1000 for ( iPhi = 0; iPhi < iPhiMax; ++iPhi )
1006 phiMax = 72*iPhi + 72;
1008 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1010 if ( firstChamber == lastChamber ) {
1011 sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d",phiMin,phiMax,firstChamber);
1012 sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d",phiMin,phiMax,lastChamber);
1015 sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
1016 sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
1019 hDeltaY [iPhi] -> SetTitle(hYtitle);
1020 hDeltaX [iPhi] -> SetTitle(hXtitle);
1022 hDeltaY [iPhi] -> Fit("fY","R","E");
1023 sigmaY [iPhi] = fY -> GetParameter(2) ;
1024 errSY [iPhi] = fY -> GetParError(2) ;
1026 rmsX [iPhi] = hDeltaX [iPhi] -> GetRMS();
1027 errSX [iPhi] = hDeltaX [iPhi] -> GetRMSError();
1029 phi [iPhi] = 72*iPhi + 36 ;
1033 //--------------------------------------------------------------------------------------------------------------
1034 //For plotting resolution in function of the angle of incidence
1036 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1040 TGraphErrors * GraphX = new TGraphErrors( iPhiMax, phi, rmsX, errPh, errSX);
1041 GraphX->SetTitle("Resolution en X (Phi)");
1042 GraphX->Draw("ALP");
1045 TGraphErrors * GraphY = new TGraphErrors( iPhiMax, phi, sigmaY, errPh, errSY);
1046 GraphY->SetTitle("Resolution en Y (Phi)");
1047 GraphY->Draw("ALP");
1051 MUONLoader->UnloadTracks();
1061 /*****************************************************************************************************************/
1062 /*****************************************************************************************************************/
1063 /*RESOLUTION IN FUNCTION OF THETA*/
1065 void resolutionTheta( Int_t event2Check=0, char * filename="galice.root" )
1067 cout<<"\nChamber(s) chosen;\nFirst chamber:";
1069 cout<<"Last chamber:";
1076 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
1077 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
1078 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
1079 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
1080 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
1081 hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
1082 hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
1083 hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
1084 hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
1086 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
1087 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
1088 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
1089 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
1090 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
1091 hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
1092 hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
1093 hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
1094 hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
1097 //Creating Run Loader and openning file containing Hits
1098 //--------------------------------------------------------------------------------------------------------------
1099 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
1100 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
1103 //Getting MUONLoader
1104 //--------------------------------------------------------------------------------------------------------------
1105 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1106 MUONLoader -> LoadTracks("READ");
1107 MUONLoader -> LoadRecPoints("READ");
1110 //Creating a MUON data container
1111 //--------------------------------------------------------------------------------------------------------------
1112 AliMUONData muondata(MUONLoader,"MUON","MUON");
1115 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1118 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1120 if ( event2Check!=0 )
1122 printf("\r>>> Event %d ",iEvent);
1123 RunLoader->GetEvent(iEvent);
1126 muondata.SetTreeAddress("RT");
1127 muondata.GetRecTracks();
1128 tracks = muondata.RecTracks();
1131 nTracks = (Int_t) tracks -> GetEntriesFast();
1132 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1134 track = (AliMUONTrack*) tracks ->At(iTrack) ;
1135 trackParams = track ->GetTrackParamAtCluster();
1136 hits = track ->GetHitForRecAtHit() ;
1139 nHits = (Int_t) hits -> GetEntriesFast();
1140 for ( iHit = 0; iHit < nHits; ++iHit )
1142 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
1143 hit = (AliMUONHitForRec* ) hits -> At(iHit);
1144 chamberNbr = hit -> GetChamberNumber();
1145 if ( chamberNbr >= firstChamber -1 ) {
1146 if ( chamberNbr < lastChamber ) {
1148 Double_t traX, traY;
1149 Double_t hitX, hitY, hitZ;
1150 Double_t aveY, aveX;
1160 traX = trackParam -> GetNonBendingCoor();
1161 traY = trackParam -> GetBendingCoor() ;
1163 hitX = hit -> GetNonBendingCoor();
1164 hitY = hit -> GetBendingCoor() ;
1165 hitZ = hit -> GetZ();
1167 aveX = 1./2.*(traX + hitX);
1168 aveY = 1./2.*(traY + hitY);
1170 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
1171 theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ ));
1173 iTheta = int( theta );
1174 if( theta > 10 ) iTheta = 9;
1175 if( theta < 1 ) iTheta = 1;
1177 deltaX = traX - hitX;
1178 deltaY = traY - hitY;
1180 hDeltaX [iTheta-1] -> Fill( deltaX );
1181 hDeltaY [iTheta-1] -> Fill( deltaY );
1189 //End loop on tracks
1190 muondata.ResetRecTracks();
1191 if ( event2Check!=0 )
1195 //End loop on events
1199 Int_t iThetaMax = 9;
1200 Int_t thetaMin, thetaMax;
1209 for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1214 //To find the polar angle
1215 thetaMin = 178 - iTheta;
1216 thetaMax = 179 - iTheta;
1218 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1220 if ( firstChamber == lastChamber ) {
1221 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
1222 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
1225 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1226 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1229 hDeltaY [iTheta] -> SetTitle(hYtitle);
1230 hDeltaX [iTheta] -> SetTitle(hXtitle);
1232 hDeltaY [iTheta] -> Fit("fY","R","E");
1233 sigmaY [iTheta] = fY -> GetParameter(2) ;
1234 errSY [iTheta] = fY -> GetParError(2) ;
1236 rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS();
1237 errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1239 theta [iTheta] = 178.5 - iTheta;
1240 ErrTh [iTheta] = 0.5;
1243 //--------------------------------------------------------------------------------------------------------------
1244 //For plotting resolution in function of the angle of incidence
1246 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1250 TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, ErrTh, errSX);
1251 GraphX->SetTitle("Resolution en X (Theta)");
1252 GraphX->Draw("ALP");
1255 TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, ErrTh, errSY);
1256 GraphY->SetTitle("Resolution en Y (Theta)");
1257 GraphY->Draw("ALP");
1260 MUONLoader->UnloadTracks();
1268 /*****************************************************************************************************************/
1269 /*****************************************************************************************************************/
1270 /*RESOLUTION IN FUNCTION OF THETA OF INCIDENCE*/
1272 void resolutionThetaI( Int_t event2Check=0, char * filename="galice.root" )
1274 cout<<"\nChamber(s) chosen;\nFirst chamber:";
1276 cout<<"Last chamber:";
1283 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
1284 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
1285 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
1286 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
1287 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
1288 hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
1289 hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
1290 hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
1291 hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
1292 hDeltaX[9] = new TH1F ("hDeltaX9" , "No Interest", 100, -0.3, 0.3);
1293 hDeltaX[10]= new TH1F ("hDeltaX10", "No Interest", 100, -0.3, 0.3);
1295 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
1296 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
1297 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
1298 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
1299 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
1300 hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
1301 hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
1302 hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
1303 hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
1304 hDeltaY[9] = new TH1F ("hDeltaY9" , "No Interest", 500, -0.1, 0.1);
1305 hDeltaY[10]= new TH1F ("hDeltaY10", "No Interest", 500, -0.1, 0.1);
1308 //Creating Run Loader and openning file containing Hits
1309 //--------------------------------------------------------------------------------------------------------------
1310 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
1311 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
1314 //Getting MUONLoader
1315 //--------------------------------------------------------------------------------------------------------------
1316 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1317 MUONLoader -> LoadTracks("READ");
1318 MUONLoader -> LoadRecPoints("READ");
1321 //Creating a MUON data container
1322 //--------------------------------------------------------------------------------------------------------------
1323 AliMUONData muondata(MUONLoader,"MUON","MUON");
1326 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1329 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1331 if ( event2Check!=0 )
1333 printf("\r>>> Event %d ",iEvent);
1334 RunLoader->GetEvent(iEvent);
1337 muondata.SetTreeAddress("RT");
1338 muondata.GetRecTracks();
1339 tracks = muondata.RecTracks();
1342 nTracks = (Int_t) tracks -> GetEntriesFast();
1343 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1345 track = (AliMUONTrack*) tracks ->At(iTrack) ;
1346 trackParams = track ->GetTrackParamAtCluster();
1347 hits = track ->GetHitForRecAtHit() ;
1350 nHits = (Int_t) hits -> GetEntriesFast();
1351 for ( iHit = 0; iHit < nHits; ++iHit )
1353 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
1354 hit = (AliMUONHitForRec* ) hits -> At(iHit);
1355 chamberNbr = hit -> GetChamberNumber();
1356 if ( chamberNbr >= firstChamber - 1 ) {
1357 if ( chamberNbr < lastChamber ) {
1359 Double_t Px, Py, Pz, Pr;
1362 Double_t traX, traY;
1363 Double_t hitX, hitY;
1373 Px = trackParam -> Px() ;
1374 Py = trackParam -> Py() ;
1375 Pz = trackParam -> Pz() ;
1376 Pr = sqrt( Px*Px + Py*Py );
1378 traX = trackParam -> GetNonBendingCoor();
1379 traY = trackParam -> GetBendingCoor();
1381 hitX = hit -> GetNonBendingCoor();
1382 hitY = hit -> GetBendingCoor();
1384 //The calculation of theta, the angle of incidence:
1385 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr))
1387 if ( theta < 79 ) iTheta = 0;
1388 else iTheta = int( theta - 79.);
1390 deltaX = traX - hitX;
1391 deltaY = traY - hitY;
1393 hDeltaX [iTheta] -> Fill (deltaX);
1394 hDeltaY [iTheta] -> Fill (deltaY);
1402 //End loop on tracks
1403 muondata.ResetRecTracks();
1404 if ( event2Check!=0 )
1408 //End loop on events
1412 Int_t iThetaMax = 11;
1413 Int_t thetaMin, thetaMax;
1422 for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1427 thetaMin = iTheta + 79;
1428 thetaMax = iTheta + 80;
1430 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1432 if ( firstChamber == lastChamber ) {
1433 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
1434 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
1437 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1438 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1441 hDeltaY [iTheta] -> SetTitle(hYtitle);
1442 hDeltaX [iTheta] -> SetTitle(hXtitle);
1444 hDeltaY [iTheta] -> Fit("fY","R","E");
1445 sigmaY [iTheta] = fY -> GetParameter(2) ;
1446 errSY [iTheta] = fY -> GetParError(2) ;
1448 rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS();
1449 errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1451 theta [iTheta] = iTheta + 79.5 ;
1452 errTh [iTheta] = 0.5;
1455 //--------------------------------------------------------------------------------------------------------------
1456 //For plotting resolution in function of the angle of incidence
1458 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1462 TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, errTh, errSX);
1463 GraphX->SetTitle("Resolution en X (Theta)");
1464 GraphX->Draw("ALP");
1467 TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, errTh, errSY);
1468 GraphY->SetTitle("Resolution en Y (Theta)");
1469 GraphY->Draw("ALP");
1472 MUONLoader->UnloadTracks();