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 **************************************************************************/
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
26 #if !defined(__CINT__) || defined(__MAKECINT__)
30 #include "TClonesArray.h"
33 #include "TParticle.h"
44 #include "TGraphErrors.h"
52 #include "AliRunLoader.h"
53 #include "AliHeader.h"
54 #include "AliLoader.h"
55 #include "AliTracker.h"
57 #include "AliMagFMaps.h"
62 #include "AliMUONData.h"
63 #include "AliMUONConstants.h"
65 #include "AliMUONHit.h"
66 #include "AliMUONHitForRec.h"
67 #include "AliMUONTrack.h"
68 #include "AliMUONTrackParam.h"
69 #include "AliMUONTrackExtrap.h"
71 #include "AliMpVSegmentation.h"
72 #include "AliMpIntPair.h"
73 #include "AliMpDEManager.h"
78 const Double_t kInvPi = 1./3.14159265;
84 Int_t nEvents, iEvent;
86 Int_t nTracks, iTrack;
90 Int_t firstChamber, lastChamber;
93 AliMUONTrack * track ;
94 AliMUONHitForRec * hit = 0;
95 AliMUONTrackParam * trackParam = 0;
97 TClonesArray * tracks ;
98 TClonesArray * trackParams;
106 /*****************************************************************************************************************/
107 /*****************************************************************************************************************/
110 void efficiency( Int_t event2Check=0, char * filename="galice.root" )
112 cout<<"\nChamber(s) chosen;\nFirst chamber:";
114 cout<<"Last chamber:";
118 //Creating Run Loader and openning file containing Hits
119 //--------------------------------------------------------------------------------------------------------------
120 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
121 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
125 //--------------------------------------------------------------------------------------------------------------
126 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
127 MUONLoader -> LoadTracks("READ");
128 MUONLoader -> LoadRecPoints("READ");
131 //Creating a MUON data container
132 //--------------------------------------------------------------------------------------------------------------
133 AliMUONData muondata(MUONLoader,"MUON","MUON");
136 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
140 Int_t chamber[10] = {0};
141 Int_t detEltNew, detElt;
144 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
146 if ( event2Check!=0 )
148 printf("\r>>> Event %d ",iEvent);
149 RunLoader -> GetEvent(iEvent);
151 muondata.SetTreeAddress("RT");
152 muondata.GetRecTracks();
153 tracks = muondata.RecTracks();
156 nTracks = (Int_t) tracks -> GetEntriesFast();
157 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
159 track = (AliMUONTrack*) tracks -> At(iTrack);
160 hits = track -> GetHitForRecAtHit();
163 nHits = (Int_t) hits -> GetEntriesFast();
165 for ( iHit = 0; iHit < nHits; ++iHit )
167 hit = (AliMUONHitForRec*) hits -> At(iHit);
168 chamberNbr = hit -> GetChamberNumber();
169 detElt = hit -> GetDetElemId();
170 detEltNew = int(detElt/100);
171 if( chamberNbr >= firstChamber-1 ) {
172 if( chamberNbr < lastChamber ) {
173 if( detEltNew == detEltOld )
176 chamber[chamberNbr] = chamber[chamberNbr] + 1;
177 detEltOld = detEltNew;
186 muondata.ResetRecTracks();
187 if (event2Check != 0)
189 trackNb = trackNb + nTracks;
192 //--------------------------------------------------------------------------------------------------------------
195 for (Int_t i = firstChamber-1; i < lastChamber; ++i )
197 printf ("\nChamber %d has responded %d times on %d tracks", i+1, chamber[i], trackNb);
199 cout<<"\nEfficiency = ? (IS UNKNOWN)\n";
201 Double_t eff = chamber[i]*100./trackNb; cout<<"\nEfficiency = "<<eff<<" %\n";
205 MUONLoader->UnloadTracks();
212 /*****************************************************************************************************************/
213 /*****************************************************************************************************************/
214 /*EFFICIENCY IN FUNCTION OF THETA AND PHI*/
216 void efficiencyThetaPhi( Int_t event2Check=0, char * filename="galice.root" )
218 cout<<"\nChamber(s) chosen;\nFirst chamber:";
220 cout<<"Last chamber:";
224 Int_t eff [54] = {0};
225 Int_t trackN [54] = {0};
229 //Creating Run Loader and openning file containing Hits
230 //--------------------------------------------------------------------------------------------------------------
231 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
232 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
236 //--------------------------------------------------------------------------------------------------------------
237 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
238 MUONLoader -> LoadTracks("READ");
239 MUONLoader -> LoadRecPoints("READ");
242 //--------------------------------------------------------------------------------------------------------------
243 //Set mag field; waiting for mag field in CDB
244 printf("Loading field map...\n");
245 if (!AliTracker::GetFieldMap()) {
246 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
247 AliTracker::SetFieldMap(field, kFALSE);}
250 //--------------------------------------------------------------------------------------------------------------
251 //Set Field Map for track extrapolation
252 AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
255 //Creating a MUON data container
256 //--------------------------------------------------------------------------------------------------------------
257 AliMUONData muondata(MUONLoader,"MUON","MUON");
260 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
263 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
265 if ( event2Check!=0 )
267 printf("\r>>> Event %d ",iEvent);
268 RunLoader->GetEvent(iEvent);
271 muondata.SetTreeAddress("RT");
272 muondata.GetRecTracks();
273 tracks = muondata.RecTracks();
276 nTracks = (Int_t) tracks -> GetEntriesFast();
277 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
279 track = (AliMUONTrack*) tracks ->At(iTrack) ;
280 trackParams = track ->GetTrackParamAtCluster();
281 hits = track ->GetHitForRecAtHit() ;
282 chamber = firstChamber-1;
287 nHits = (Int_t) hits->GetEntriesFast();
288 for ( iHit = 0; iHit<nHits; ++iHit )
290 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
291 hit = (AliMUONHitForRec* ) hits -> At(iHit);
292 chamberNbr = hit -> GetChamberNumber();
294 if ( chamberNbr == oldChamber )
297 oldChamber = chamberNbr;
298 if ( chamberNbr > chamber - k ) {
299 if ( chamber < lastChamber ) {
300 if ( chamberNbr == chamber ) {
302 Double_t traX, traY, traZ;
303 Double_t hitX, hitY, hitZ;
304 Double_t aveX, aveY ;
312 traX = trackParam -> GetNonBendingCoor();
313 traY = trackParam -> GetBendingCoor() ;
314 traZ = trackParam -> GetZ() ;
316 hitX = hit -> GetNonBendingCoor();
317 hitY = hit -> GetBendingCoor() ;
318 hitZ = hit -> GetZ() ;
320 aveX = 1./2.*(traX + hitX);
321 aveY = 1./2.*(traY + hitY);
323 //The calculation of phi:
324 phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));
326 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
327 theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ ));
329 if ( phi < 0.) phi = 360 - fabs(phi);
330 iPhi = int( phi/72. );
331 iTheta = int( theta );
332 if( theta > 10 ) iTheta = 9;
333 if( theta < 1 ) iTheta = 1;
335 eff [9*iPhi+iTheta-1] = eff [9*iPhi+iTheta-1] + 1;
336 trackN [9*iPhi+iTheta-1] = trackN [9*iPhi+iTheta-1] + 1;
337 chamber = chamber + 1;
342 Double_t chamberZpos;
343 Double_t exXpos, exYpos;
351 chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
352 AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
353 exXpos = (Double_t) trackParam->GetNonBendingCoor();
354 exYpos = (Double_t) trackParam->GetBendingCoor();
356 //The calculation of phi:
357 phi = 180. * kInvPi * (TMath::ATan2( exYpos, exXpos ));
359 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
360 theta = 180. * kInvPi * (TMath::ATan2( sqrt(exXpos*exXpos+exYpos*exYpos), - chamberZpos ));
362 if ( phi < 0.) phi = 360. - fabs(phi);
363 iPhi = int( phi/72. );
364 iTheta = int( theta );
365 if( theta > 10 ) iTheta = 9;
366 if( theta < 1 ) iTheta = 1;
368 eff [9*iPhi+iTheta-1] = eff [9*iPhi+iTheta-1] + 0;
369 trackN [9*iPhi+iTheta-1] = trackN [9*iPhi+iTheta-1] + 1;
370 chamber = chamber + 1;
377 if ( iHit == nHits-1 ) {
378 if ( chamber == lastChamber )
392 muondata.ResetRecTracks();
393 if ( event2Check!=0 )
399 TCanvas * c1 = new TCanvas();
400 TGraph2D* effPhiTheta = new TGraph2D();
401 Double_t efficiency = 0;
403 if ( firstChamber == lastChamber )
405 for ( Int_t ph = 0; ph < 5; ++ph )
407 for ( Int_t th = 1; th < 10; ++th )
410 cout<<"\nFor Phi = ["<<ph*72<<","<<ph*72+72<<"] and Theta = ["<<179-th<<","<<180-th<<"]:";
411 cout<<"\nThe chamber "<<firstChamber<<" has responded "<<eff[i]<<" times on "<<trackN[i]<<" tracks\n";
413 Double_t e = eff [i] ;
414 Double_t n = trackN[i] ;
415 Double_t p = ph*72.+36.;
416 Double_t t = th*1. +0.5;
418 if ( trackN[i] == 0 ) {
420 cout<<"Efficiency = ? % ( IS UNKNOWN )\n";
423 efficiency = 100.*e/n;
424 cout<<"Efficiency = "<<efficiency<<" %\n";
427 effPhiTheta -> SetPoint( i, p, t, efficiency);
433 for ( Int_t ph = 0; ph < 5; ++ph )
435 for ( Int_t th = 1; th < 10; ++th )
438 cout<<"\nFor Phi = ["<<ph*72<<","<<ph*72+72<<"] and Theta = ["<<179-th<<","<<180-th<<"]:";
439 cout<<"\nThe chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff[i]<<" times on "<<trackN[i]<<" tracks\n";
441 Double_t e = eff [i] ;
442 Double_t n = trackN[i] ;
443 Double_t p = ph*72.+36.;
444 Double_t t = th*1. +0.5;
446 if ( trackN[i] == 0 ) {
448 cout<<"Efficiency = ? % ( IS UNKNOWN )\n";
451 efficiency = 100.*e/n;
452 cout<<"Efficiency = "<<efficiency<<" %\n";
455 effPhiTheta -> SetPoint( i, p, t, efficiency);
460 gStyle->SetPalette(1);
461 effPhiTheta -> Draw("surf1");
464 MUONLoader->UnloadTracks();
475 /*****************************************************************************************************************/
476 /*****************************************************************************************************************/
477 /*EFFICIENCY IN FUNCTION OF THETA OF INCIDENCE*/
480 void efficiencyThetaI( Int_t event2Check=0, char * filename="galice.root" )
482 cout<<"\nChamber(s) chosen;\nFirst chamber:";
484 cout<<"Last chamber:";
488 Int_t eff [12] = {0};
489 Int_t trackN [12] = {0};
494 //Creating Run Loader and openning file containing Hits
495 //--------------------------------------------------------------------------------------------------------------
496 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
497 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
501 //--------------------------------------------------------------------------------------------------------------
502 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
503 MUONLoader -> LoadTracks("READ");
504 MUONLoader -> LoadRecPoints("READ");
507 //--------------------------------------------------------------------------------------------------------------
508 //Set mag field; waiting for mag field in CDB
509 printf("Loading field map...\n");
510 if (!AliTracker::GetFieldMap()) {
511 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
512 AliTracker::SetFieldMap(field, kFALSE);}
515 //--------------------------------------------------------------------------------------------------------------
516 //Set Field Map for track extrapolation
517 AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
520 //Creating a MUON data container
521 //--------------------------------------------------------------------------------------------------------------
522 AliMUONData muondata(MUONLoader,"MUON","MUON");
525 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
528 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
530 if ( event2Check!=0 )
532 printf("\r>>> Event %d ",iEvent);
533 RunLoader->GetEvent(iEvent);
536 muondata.SetTreeAddress("RT");
537 muondata.GetRecTracks();
538 tracks = muondata.RecTracks();
541 nTracks = (Int_t) tracks -> GetEntriesFast();
542 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
544 track = (AliMUONTrack*) tracks ->At(iTrack) ;
545 trackParams = track ->GetTrackParamAtCluster();
546 hits = track ->GetHitForRecAtHit() ;
547 chamber = firstChamber - 1;
552 nHits = (Int_t) hits -> GetEntriesFast();
553 for ( iHit = 0; iHit < nHits; ++iHit )
555 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
556 hit = (AliMUONHitForRec*) hits -> At(iHit);
557 chamberNbr = hit -> GetChamberNumber();
559 if ( chamberNbr == oldChamber )
562 oldChamber = chamberNbr;
563 if ( chamberNbr > chamber - k ) {
564 if ( chamber < lastChamber ) {
565 if ( chamberNbr == chamber ) {
567 Double_t Px, Py, Pz, Pr;
573 Px = trackParam -> Px() ;
574 Py = trackParam -> Py() ;
575 Pz = trackParam -> Pz() ;
576 Pr = sqrt( Px*Px + Py*Py );
578 //The calculation of theta, the angle of incidence:
579 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
581 if ( theta < 79 ) iTheta = 11;
582 else if ( theta < 90 ) iTheta = int( theta - 79.);
585 eff [iTheta] = eff [iTheta] + 1;
586 trackN [iTheta] = trackN [iTheta] + 1;
587 chamber = chamber + 1;
592 Double_t chamberZpos;
595 Double_t Px, Py, Pz, Pr;
601 chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
602 AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
604 Px = trackParam -> Px() ;
605 Py = trackParam -> Py() ;
606 Pz = trackParam -> Pz() ;
607 Pr = sqrt( Px*Px + Py*Py );
609 //The calculation of thetaI, the angle of incidence:
610 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
612 if ( theta < 79 ) iTheta = 11;
613 else if ( theta < 90 ) iTheta = int( theta - 79.);
616 eff [iTheta] = eff [iTheta] + 0;
617 trackN [iTheta] = trackN [iTheta] + 1;
618 chamber = chamber + 1;
625 if ( iHit == nHits-1 ) {
626 if ( chamber == lastChamber )
640 muondata.ResetRecTracks();
641 if ( event2Check!=0 )
647 Double_t efficiency [11];
651 TGraph * effTheta = new TGraph ();
653 if ( firstChamber == lastChamber ) {
654 for ( th = 0; th < 11; ++th )
656 cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
657 cout<<"The chamber "<<firstChamber<<" has responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
660 Double_t e = eff [th];
661 Double_t n = trackN [th];
664 efficiency [th] = 0.;
665 cout<<"Efficiency = ? % (IS UNKNOWN) \n";
668 efficiency [th] = 100.*e/n;
669 cout<<"Efficiency = "<<efficiency [th]<<" %\n";
675 for ( th = 0; th < 11; ++th )
677 cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
678 cout<<"The chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
681 Double_t e = eff [th];
682 Double_t n = trackN [th];
685 efficiency [th] = 0.;
686 cout<<"Efficiency = ? % (IS UNKNOWN) \n";
689 efficiency [th] = 100.*e/n;
690 cout<<"Efficiency = "<<efficiency [th]<<" %\n";
695 TCanvas * c1 = new TCanvas ();
696 effTheta = new TGraph( i, t, efficiency );
698 effTheta -> Draw("ALP");
700 MUONLoader->UnloadTracks();
709 /*****************************************************************************************************************/
710 /*****************************************************************************************************************/
713 void resolution( Int_t event2Check=0, char * filename="galice.root" )
715 cout<<"\nChamber(s) chosen;\nFirst chamber:";
717 cout<<"Last chamber:";
725 hDeltax = new TH1F ("hDeltax " , "No Interest", 100, -0.3, 0.3);
726 hDeltay = new TH1F ("hDeltay " , "No Interest", 500, -0.1, 0.1);
727 hDelta3D = new TH2F ("hDelta3D" , "No Interest", 100, -0.3, 0.3, 500, -0.1, 0.1);
730 //Creating Run Loader and openning file containing Hits
731 //--------------------------------------------------------------------------------------------------------------
732 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
733 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
737 //--------------------------------------------------------------------------------------------------------------
738 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
739 MUONLoader -> LoadTracks("READ");
740 MUONLoader -> LoadRecPoints("READ");
743 //Creating a MUON data container
744 //--------------------------------------------------------------------------------------------------------------
745 AliMUONData muondata(MUONLoader,"MUON","MUON");
748 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
751 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
755 printf("\r>>> Event %d ",iEvent);
756 RunLoader->GetEvent(iEvent);
759 muondata.SetTreeAddress("RT");
760 muondata.GetRecTracks();
761 tracks = muondata.RecTracks();
764 nTracks = (Int_t) tracks -> GetEntriesFast();
765 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
767 track = (AliMUONTrack*) tracks ->At(iTrack) ;
768 trackParams = track ->GetTrackParamAtCluster();
769 hits = track ->GetHitForRecAtHit() ;
772 nHits = (Int_t) hits -> GetEntriesFast();
773 for ( iHit = 0; iHit < nHits; ++iHit )
775 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
776 hit = (AliMUONHitForRec*) hits -> At(iHit);
777 chamberNbr = hit -> GetChamberNumber();
778 if ( chamberNbr >= firstChamber-1 ) {
779 if ( chamberNbr < lastChamber ) {
785 Double_t deltaX, deltaY;
787 traX = trackParam -> GetNonBendingCoor();
788 traY = trackParam -> GetBendingCoor();
789 hitX = hit -> GetNonBendingCoor();
790 hitY = hit -> GetBendingCoor();
792 deltaX = traX - hitX;
793 deltaY = traY - hitY;
795 hDeltax -> Fill (deltaX);
796 hDeltay -> Fill (deltaY);
797 hDelta3D -> Fill (deltaX, deltaY);
804 muondata.ResetRecTracks();
805 if ( event2Check!=0 )
814 if ( firstChamber == lastChamber ) {
815 sprintf(hXtitle,"DeltaX Chamber %d",firstChamber);
816 sprintf(hYtitle,"DeltaY Chamber %d",lastChamber);
817 sprintf(h3title,"DeltaX and Delta Y Chamber %d",lastChamber);
820 sprintf(hXtitle,"DeltaX Chamber %d to %d",firstChamber,lastChamber);
821 sprintf(hYtitle,"DeltaY Chamber %d to %d",firstChamber,lastChamber);
822 sprintf(h3title,"DeltaX and Delta Y Chamber %d to %d",firstChamber,lastChamber);
826 hDeltax -> SetTitle (hXtitle);
827 hDeltay -> SetTitle (hYtitle);
828 hDelta3D -> SetTitle (h3title);
830 Double_t rmsX = hDeltax -> GetRMS ();
831 Double_t errX = hDeltax -> GetRMSError();
833 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
834 hDeltay -> Fit("fY","R","E");
835 Double_t sigY = fY -> GetParameter(2);
836 Double_t errY = fY -> GetParError (2);
838 if ( firstChamber == lastChamber ) {
839 cout<<"\nThe resolution (Non Bending plane) of the chamber "<<firstChamber<<" = "<<rmsX<<" cm +-"<<errX;
840 cout<<"\nThe resolution (Bending plane) of the chamber "<<firstChamber<<" = "<<sigY<<" cm +- "<<errY;
843 cout<<"\nThe resolution (Non Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<rmsX<<" cm +- "<<errX;
844 cout<<"\nThe resolution (Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<sigY<<" cm +- "<<errY;
848 MUONLoader->UnloadTracks();
858 /*****************************************************************************************************************/
859 /*****************************************************************************************************************/
860 /*RESOLUTION IN FUNCTION OF PHI*/
862 void resolutionPhi( Int_t event2Check=0, char * filename="galice.root" )
864 cout<<"\nChamber(s) chosen;\nFirst chamber:";
866 cout<<"Last chamber:";
873 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
874 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
875 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
876 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
877 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
879 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
880 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
881 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
882 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
883 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
886 //Creating Run Loader and openning file containing Hits
887 //--------------------------------------------------------------------------------------------------------------
888 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
889 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
893 //--------------------------------------------------------------------------------------------------------------
894 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
895 MUONLoader -> LoadTracks("READ");
896 MUONLoader -> LoadRecPoints("READ");
899 //Creating a MUON data container
900 //--------------------------------------------------------------------------------------------------------------
901 AliMUONData muondata(MUONLoader,"MUON","MUON");
904 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
907 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
909 if ( event2Check!=0 )
911 printf("\r>>> Event %d ",iEvent);
912 RunLoader->GetEvent(iEvent);
915 muondata.SetTreeAddress("RT");
916 muondata.GetRecTracks();
917 tracks = muondata.RecTracks();
920 nTracks = (Int_t) tracks -> GetEntriesFast();
921 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
923 track = (AliMUONTrack*) tracks ->At(iTrack) ;
924 trackParams = track ->GetTrackParamAtCluster();
925 hits = track ->GetHitForRecAtHit() ;
928 nHits = (Int_t) hits -> GetEntriesFast();
929 for ( iHit = 0; iHit < nHits; ++iHit )
931 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
932 hit = (AliMUONHitForRec* ) hits -> At(iHit);
933 chamberNbr = hit -> GetChamberNumber();
934 if ( chamberNbr >= firstChamber -1 ) {
935 if ( chamberNbr < lastChamber ) {
949 traX = trackParam -> GetNonBendingCoor();
950 traY = trackParam -> GetBendingCoor() ;
952 hitX = hit -> GetNonBendingCoor();
953 hitY = hit -> GetBendingCoor() ;
955 aveX = 1./2.*(traX + hitX);
956 aveY = 1./2.*(traY + hitY);
958 //The calculation of phi:
959 phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));
961 if ( phi < 0.) phi = 360. - fabs(phi);
962 iPhi = int( phi/72. );
964 deltaX = traX - hitX;
965 deltaY = traY - hitY;
967 hDeltaX [iPhi] -> Fill( deltaX );
968 hDeltaY [iPhi] -> Fill( deltaY );
977 muondata.ResetRecTracks();
978 if ( event2Check!=0 )
987 Int_t phiMin, phiMax;
996 for ( iPhi = 0; iPhi < iPhiMax; ++iPhi )
1002 phiMax = 72*iPhi + 72;
1004 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1006 if ( firstChamber == lastChamber ) {
1007 sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d",phiMin,phiMax,firstChamber);
1008 sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d",phiMin,phiMax,lastChamber);
1011 sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
1012 sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
1015 hDeltaY [iPhi] -> SetTitle(hYtitle);
1016 hDeltaX [iPhi] -> SetTitle(hXtitle);
1018 hDeltaY [iPhi] -> Fit("fY","R","E");
1019 sigmaY [iPhi] = fY -> GetParameter(2) ;
1020 errSY [iPhi] = fY -> GetParError(2) ;
1022 rmsX [iPhi] = hDeltaX [iPhi] -> GetRMS();
1023 errSX [iPhi] = hDeltaX [iPhi] -> GetRMSError();
1025 phi [iPhi] = 72*iPhi + 36 ;
1029 //--------------------------------------------------------------------------------------------------------------
1030 //For plotting resolution in function of the angle of incidence
1032 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1036 TGraphErrors * GraphX = new TGraphErrors( iPhiMax, phi, rmsX, errPh, errSX);
1037 GraphX->SetTitle("Resolution en X (Phi)");
1038 GraphX->Draw("ALP");
1041 TGraphErrors * GraphY = new TGraphErrors( iPhiMax, phi, sigmaY, errPh, errSY);
1042 GraphY->SetTitle("Resolution en Y (Phi)");
1043 GraphY->Draw("ALP");
1047 MUONLoader->UnloadTracks();
1057 /*****************************************************************************************************************/
1058 /*****************************************************************************************************************/
1059 /*RESOLUTION IN FUNCTION OF THETA*/
1061 void resolutionTheta( Int_t event2Check=0, char * filename="galice.root" )
1063 cout<<"\nChamber(s) chosen;\nFirst chamber:";
1065 cout<<"Last chamber:";
1072 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
1073 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
1074 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
1075 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
1076 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
1077 hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
1078 hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
1079 hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
1080 hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
1082 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
1083 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
1084 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
1085 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
1086 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
1087 hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
1088 hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
1089 hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
1090 hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
1093 //Creating Run Loader and openning file containing Hits
1094 //--------------------------------------------------------------------------------------------------------------
1095 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
1096 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
1099 //Getting MUONLoader
1100 //--------------------------------------------------------------------------------------------------------------
1101 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1102 MUONLoader -> LoadTracks("READ");
1103 MUONLoader -> LoadRecPoints("READ");
1106 //Creating a MUON data container
1107 //--------------------------------------------------------------------------------------------------------------
1108 AliMUONData muondata(MUONLoader,"MUON","MUON");
1111 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1114 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1116 if ( event2Check!=0 )
1118 printf("\r>>> Event %d ",iEvent);
1119 RunLoader->GetEvent(iEvent);
1122 muondata.SetTreeAddress("RT");
1123 muondata.GetRecTracks();
1124 tracks = muondata.RecTracks();
1127 nTracks = (Int_t) tracks -> GetEntriesFast();
1128 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1130 track = (AliMUONTrack*) tracks ->At(iTrack) ;
1131 trackParams = track ->GetTrackParamAtCluster();
1132 hits = track ->GetHitForRecAtHit() ;
1135 nHits = (Int_t) hits -> GetEntriesFast();
1136 for ( iHit = 0; iHit < nHits; ++iHit )
1138 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
1139 hit = (AliMUONHitForRec* ) hits -> At(iHit);
1140 chamberNbr = hit -> GetChamberNumber();
1141 if ( chamberNbr >= firstChamber -1 ) {
1142 if ( chamberNbr < lastChamber ) {
1144 Double_t traX, traY;
1145 Double_t hitX, hitY, hitZ;
1146 Double_t aveY, aveX;
1156 traX = trackParam -> GetNonBendingCoor();
1157 traY = trackParam -> GetBendingCoor() ;
1159 hitX = hit -> GetNonBendingCoor();
1160 hitY = hit -> GetBendingCoor() ;
1161 hitZ = hit -> GetZ();
1163 aveX = 1./2.*(traX + hitX);
1164 aveY = 1./2.*(traY + hitY);
1166 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
1167 theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ ));
1169 iTheta = int( theta );
1170 if( theta > 10 ) iTheta = 9;
1171 if( theta < 1 ) iTheta = 1;
1173 deltaX = traX - hitX;
1174 deltaY = traY - hitY;
1176 hDeltaX [iTheta-1] -> Fill( deltaX );
1177 hDeltaY [iTheta-1] -> Fill( deltaY );
1185 //End loop on tracks
1186 muondata.ResetRecTracks();
1187 if ( event2Check!=0 )
1191 //End loop on events
1195 Int_t iThetaMax = 9;
1196 Int_t thetaMin, thetaMax;
1205 for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1210 //To find the polar angle
1211 thetaMin = 178 - iTheta;
1212 thetaMax = 179 - iTheta;
1214 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1216 if ( firstChamber == lastChamber ) {
1217 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
1218 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
1221 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1222 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1225 hDeltaY [iTheta] -> SetTitle(hYtitle);
1226 hDeltaX [iTheta] -> SetTitle(hXtitle);
1228 hDeltaY [iTheta] -> Fit("fY","R","E");
1229 sigmaY [iTheta] = fY -> GetParameter(2) ;
1230 errSY [iTheta] = fY -> GetParError(2) ;
1232 rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS();
1233 errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1235 theta [iTheta] = 178.5 - iTheta;
1236 ErrTh [iTheta] = 0.5;
1239 //--------------------------------------------------------------------------------------------------------------
1240 //For plotting resolution in function of the angle of incidence
1242 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1246 TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, ErrTh, errSX);
1247 GraphX->SetTitle("Resolution en X (Theta)");
1248 GraphX->Draw("ALP");
1251 TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, ErrTh, errSY);
1252 GraphY->SetTitle("Resolution en Y (Theta)");
1253 GraphY->Draw("ALP");
1256 MUONLoader->UnloadTracks();
1264 /*****************************************************************************************************************/
1265 /*****************************************************************************************************************/
1266 /*RESOLUTION IN FUNCTION OF THETA OF INCIDENCE*/
1268 void resolutionThetaI( Int_t event2Check=0, char * filename="galice.root" )
1270 cout<<"\nChamber(s) chosen;\nFirst chamber:";
1272 cout<<"Last chamber:";
1279 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
1280 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
1281 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
1282 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
1283 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
1284 hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
1285 hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
1286 hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
1287 hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
1288 hDeltaX[9] = new TH1F ("hDeltaX9" , "No Interest", 100, -0.3, 0.3);
1289 hDeltaX[10]= new TH1F ("hDeltaX10", "No Interest", 100, -0.3, 0.3);
1291 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
1292 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
1293 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
1294 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
1295 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
1296 hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
1297 hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
1298 hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
1299 hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
1300 hDeltaY[9] = new TH1F ("hDeltaY9" , "No Interest", 500, -0.1, 0.1);
1301 hDeltaY[10]= new TH1F ("hDeltaY10", "No Interest", 500, -0.1, 0.1);
1304 //Creating Run Loader and openning file containing Hits
1305 //--------------------------------------------------------------------------------------------------------------
1306 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
1307 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
1310 //Getting MUONLoader
1311 //--------------------------------------------------------------------------------------------------------------
1312 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1313 MUONLoader -> LoadTracks("READ");
1314 MUONLoader -> LoadRecPoints("READ");
1317 //Creating a MUON data container
1318 //--------------------------------------------------------------------------------------------------------------
1319 AliMUONData muondata(MUONLoader,"MUON","MUON");
1322 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1325 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1327 if ( event2Check!=0 )
1329 printf("\r>>> Event %d ",iEvent);
1330 RunLoader->GetEvent(iEvent);
1333 muondata.SetTreeAddress("RT");
1334 muondata.GetRecTracks();
1335 tracks = muondata.RecTracks();
1338 nTracks = (Int_t) tracks -> GetEntriesFast();
1339 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1341 track = (AliMUONTrack*) tracks ->At(iTrack) ;
1342 trackParams = track ->GetTrackParamAtCluster();
1343 hits = track ->GetHitForRecAtHit() ;
1346 nHits = (Int_t) hits -> GetEntriesFast();
1347 for ( iHit = 0; iHit < nHits; ++iHit )
1349 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
1350 hit = (AliMUONHitForRec* ) hits -> At(iHit);
1351 chamberNbr = hit -> GetChamberNumber();
1352 if ( chamberNbr >= firstChamber - 1 ) {
1353 if ( chamberNbr < lastChamber ) {
1355 Double_t Px, Py, Pz, Pr;
1358 Double_t traX, traY;
1359 Double_t hitX, hitY;
1369 Px = trackParam -> Px() ;
1370 Py = trackParam -> Py() ;
1371 Pz = trackParam -> Pz() ;
1372 Pr = sqrt( Px*Px + Py*Py );
1374 traX = trackParam -> GetNonBendingCoor();
1375 traY = trackParam -> GetBendingCoor();
1377 hitX = hit -> GetNonBendingCoor();
1378 hitY = hit -> GetBendingCoor();
1380 //The calculation of theta, the angle of incidence:
1381 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr))
1383 if ( theta < 79 ) iTheta = 0;
1384 else iTheta = int( theta - 79.);
1386 deltaX = traX - hitX;
1387 deltaY = traY - hitY;
1389 hDeltaX [iTheta] -> Fill (deltaX);
1390 hDeltaY [iTheta] -> Fill (deltaY);
1398 //End loop on tracks
1399 muondata.ResetRecTracks();
1400 if ( event2Check!=0 )
1404 //End loop on events
1408 Int_t iThetaMax = 11;
1409 Int_t thetaMin, thetaMax;
1418 for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1423 thetaMin = iTheta + 79;
1424 thetaMax = iTheta + 80;
1426 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1428 if ( firstChamber == lastChamber ) {
1429 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
1430 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
1433 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1434 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1437 hDeltaY [iTheta] -> SetTitle(hYtitle);
1438 hDeltaX [iTheta] -> SetTitle(hXtitle);
1440 hDeltaY [iTheta] -> Fit("fY","R","E");
1441 sigmaY [iTheta] = fY -> GetParameter(2) ;
1442 errSY [iTheta] = fY -> GetParError(2) ;
1444 rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS();
1445 errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1447 theta [iTheta] = iTheta + 79.5 ;
1448 errTh [iTheta] = 0.5;
1451 //--------------------------------------------------------------------------------------------------------------
1452 //For plotting resolution in function of the angle of incidence
1454 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1458 TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, errTh, errSX);
1459 GraphX->SetTitle("Resolution en X (Theta)");
1460 GraphX->Draw("ALP");
1463 TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, errTh, errSY);
1464 GraphY->SetTitle("Resolution en Y (Theta)");
1465 GraphY->Draw("ALP");
1468 MUONLoader->UnloadTracks();