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 "AliMUONRawCluster.h"
68 #include "AliMUONTrack.h"
69 #include "AliMUONTrackParam.h"
70 #include "AliMUONTrackExtrap.h"
72 #include "AliMpVSegmentation.h"
73 #include "AliMpIntPair.h"
74 #include "AliMpDEManager.h"
79 const Double_t kInvPi = 1./3.14159265;
85 Int_t nEvents, iEvent;
87 Int_t nTracks, iTrack;
91 Int_t firstChamber, lastChamber;
94 AliMUONTrack * track ;
95 AliMUONHitForRec * hit = 0;
96 AliMUONTrackParam * trackParam = 0;
98 TClonesArray * tracks ;
99 TClonesArray * trackParams;
100 TClonesArray * hits ;
107 /*****************************************************************************************************************/
108 /*****************************************************************************************************************/
111 void efficiency( Int_t event2Check=0, char * filename="galice.root" )
113 cout<<"\nChamber(s) chosen;\nFirst chamber:";
115 cout<<"Last chamber:";
119 //Creating Run Loader and openning file containing Hits
120 //--------------------------------------------------------------------------------------------------------------
121 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
122 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
126 //--------------------------------------------------------------------------------------------------------------
127 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
128 MUONLoader -> LoadTracks("READ");
129 MUONLoader -> LoadRecPoints("READ");
132 //Creating a MUON data container
133 //--------------------------------------------------------------------------------------------------------------
134 AliMUONData muondata(MUONLoader,"MUON","MUON");
137 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
141 Int_t chamber[10] = {0};
142 Int_t detEltNew, detElt;
145 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
147 if ( event2Check!=0 )
149 printf("\r>>> Event %d ",iEvent);
150 RunLoader -> GetEvent(iEvent);
152 muondata.SetTreeAddress("RT");
153 muondata.GetRecTracks();
154 tracks = muondata.RecTracks();
157 nTracks = (Int_t) tracks -> GetEntriesFast();
158 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
160 track = (AliMUONTrack*) tracks -> At(iTrack);
161 hits = track -> GetHitForRecAtHit();
164 nHits = (Int_t) hits -> GetEntriesFast();
166 for ( iHit = 0; iHit < nHits; ++iHit )
168 hit = (AliMUONHitForRec*) hits -> At(iHit);
169 chamberNbr = hit -> GetChamberNumber();
170 detElt = hit -> GetDetElemId();
171 detEltNew = int(detElt/100);
172 if( chamberNbr >= firstChamber-1 ) {
173 if( chamberNbr < lastChamber ) {
174 if( detEltNew == detEltOld )
177 chamber[chamberNbr] = chamber[chamberNbr] + 1;
178 detEltOld = detEltNew;
187 muondata.ResetRecTracks();
188 if (event2Check != 0)
190 trackNb = trackNb + nTracks;
193 //--------------------------------------------------------------------------------------------------------------
196 for (Int_t i = firstChamber-1; i < lastChamber; ++i )
198 printf ("\nChamber %d has responded %d times on %d tracks", i+1, chamber[i], trackNb);
200 cout<<"\nEfficiency = ? (IS UNKNOWN)\n";
202 Double_t eff = chamber[i]*100./trackNb; cout<<"\nEfficiency = "<<eff<<" %\n";
206 MUONLoader->UnloadTracks();
213 /*****************************************************************************************************************/
214 /*****************************************************************************************************************/
215 /*EFFICIENCY IN FUNCTION OF THETA AND PHI*/
217 void efficiencyThetaPhi( Int_t event2Check=0, char * filename="galice.root" )
219 cout<<"\nChamber(s) chosen;\nFirst chamber:";
221 cout<<"Last chamber:";
225 Int_t eff [54] = {0};
226 Int_t trackN [54] = {0};
230 //Creating Run Loader and openning file containing Hits
231 //--------------------------------------------------------------------------------------------------------------
232 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
233 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
237 //--------------------------------------------------------------------------------------------------------------
238 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
239 MUONLoader -> LoadTracks("READ");
240 MUONLoader -> LoadRecPoints("READ");
243 //--------------------------------------------------------------------------------------------------------------
244 //Set mag field; waiting for mag field in CDB
245 printf("Loading field map...\n");
246 if (!AliTracker::GetFieldMap()) {
247 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
248 AliTracker::SetFieldMap(field, kFALSE);}
251 //--------------------------------------------------------------------------------------------------------------
252 //Set Field Map for track extrapolation
253 AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
256 //Creating a MUON data container
257 //--------------------------------------------------------------------------------------------------------------
258 AliMUONData muondata(MUONLoader,"MUON","MUON");
261 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
264 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
266 if ( event2Check!=0 )
268 printf("\r>>> Event %d ",iEvent);
269 RunLoader->GetEvent(iEvent);
272 muondata.SetTreeAddress("RT");
273 muondata.GetRecTracks();
274 tracks = muondata.RecTracks();
277 nTracks = (Int_t) tracks -> GetEntriesFast();
278 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
280 track = (AliMUONTrack*) tracks ->At(iTrack) ;
281 trackParams = track ->GetTrackParamAtHit();
282 hits = track ->GetHitForRecAtHit() ;
283 chamber = firstChamber-1;
288 nHits = (Int_t) hits->GetEntriesFast();
289 for ( iHit = 0; iHit<nHits; ++iHit )
291 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
292 hit = (AliMUONHitForRec* ) hits -> At(iHit);
293 chamberNbr = hit -> GetChamberNumber();
295 if ( chamberNbr == oldChamber )
298 oldChamber = chamberNbr;
299 if ( chamberNbr > chamber - k ) {
300 if ( chamber < lastChamber ) {
301 if ( chamberNbr == chamber ) {
303 Double_t traX, traY, traZ;
304 Double_t hitX, hitY, hitZ;
305 Double_t aveX, aveY ;
313 traX = trackParam -> GetNonBendingCoor();
314 traY = trackParam -> GetBendingCoor() ;
315 traZ = trackParam -> GetZ() ;
317 hitX = hit -> GetNonBendingCoor();
318 hitY = hit -> GetBendingCoor() ;
319 hitZ = hit -> GetZ() ;
321 aveX = 1./2.*(traX + hitX);
322 aveY = 1./2.*(traY + hitY);
324 //The calculation of phi:
325 phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));
327 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
328 theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ ));
330 if ( phi < 0.) phi = 360 - fabs(phi);
331 iPhi = int( phi/72. );
332 iTheta = int( theta );
333 if( theta > 10 ) iTheta = 9;
334 if( theta < 1 ) iTheta = 1;
336 eff [9*iPhi+iTheta-1] = eff [9*iPhi+iTheta-1] + 1;
337 trackN [9*iPhi+iTheta-1] = trackN [9*iPhi+iTheta-1] + 1;
338 chamber = chamber + 1;
343 Double_t chamberZpos;
344 Double_t exXpos, exYpos;
352 chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
353 AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
354 exXpos = (Double_t) trackParam->GetNonBendingCoor();
355 exYpos = (Double_t) trackParam->GetBendingCoor();
357 //The calculation of phi:
358 phi = 180. * kInvPi * (TMath::ATan2( exYpos, exXpos ));
360 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
361 theta = 180. * kInvPi * (TMath::ATan2( sqrt(exXpos*exXpos+exYpos*exYpos), - chamberZpos ));
363 if ( phi < 0.) phi = 360. - fabs(phi);
364 iPhi = int( phi/72. );
365 iTheta = int( theta );
366 if( theta > 10 ) iTheta = 9;
367 if( theta < 1 ) iTheta = 1;
369 eff [9*iPhi+iTheta-1] = eff [9*iPhi+iTheta-1] + 0;
370 trackN [9*iPhi+iTheta-1] = trackN [9*iPhi+iTheta-1] + 1;
371 chamber = chamber + 1;
378 if ( iHit == nHits-1 ) {
379 if ( chamber == lastChamber )
393 muondata.ResetRecTracks();
394 if ( event2Check!=0 )
400 TCanvas * c1 = new TCanvas();
401 TGraph2D* effPhiTheta = new TGraph2D();
402 Double_t efficiency = 0;
404 if ( firstChamber == lastChamber )
406 for ( Int_t ph = 0; ph < 5; ++ph )
408 for ( Int_t th = 1; th < 10; ++th )
411 cout<<"\nFor Phi = ["<<ph*72<<","<<ph*72+72<<"] and Theta = ["<<179-th<<","<<180-th<<"]:";
412 cout<<"\nThe chamber "<<firstChamber<<" has responded "<<eff[i]<<" times on "<<trackN[i]<<" tracks\n";
414 Double_t e = eff [i] ;
415 Double_t n = trackN[i] ;
416 Double_t p = ph*72.+36.;
417 Double_t t = th*1. +0.5;
419 if ( trackN[i] == 0 ) {
421 cout<<"Efficiency = ? % ( IS UNKNOWN )\n";
424 efficiency = 100.*e/n;
425 cout<<"Efficiency = "<<efficiency<<" %\n";
428 effPhiTheta -> SetPoint( i, p, t, efficiency);
434 for ( Int_t ph = 0; ph < 5; ++ph )
436 for ( Int_t th = 1; th < 10; ++th )
439 cout<<"\nFor Phi = ["<<ph*72<<","<<ph*72+72<<"] and Theta = ["<<179-th<<","<<180-th<<"]:";
440 cout<<"\nThe chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff[i]<<" times on "<<trackN[i]<<" tracks\n";
442 Double_t e = eff [i] ;
443 Double_t n = trackN[i] ;
444 Double_t p = ph*72.+36.;
445 Double_t t = th*1. +0.5;
447 if ( trackN[i] == 0 ) {
449 cout<<"Efficiency = ? % ( IS UNKNOWN )\n";
452 efficiency = 100.*e/n;
453 cout<<"Efficiency = "<<efficiency<<" %\n";
456 effPhiTheta -> SetPoint( i, p, t, efficiency);
461 gStyle->SetPalette(1);
462 effPhiTheta -> Draw("surf1");
465 MUONLoader->UnloadTracks();
476 /*****************************************************************************************************************/
477 /*****************************************************************************************************************/
478 /*EFFICIENCY IN FUNCTION OF THETA OF INCIDENCE*/
481 void efficiencyThetaI( Int_t event2Check=0, char * filename="galice.root" )
483 cout<<"\nChamber(s) chosen;\nFirst chamber:";
485 cout<<"Last chamber:";
489 Int_t eff [12] = {0};
490 Int_t trackN [12] = {0};
495 //Creating Run Loader and openning file containing Hits
496 //--------------------------------------------------------------------------------------------------------------
497 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
498 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
502 //--------------------------------------------------------------------------------------------------------------
503 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
504 MUONLoader -> LoadTracks("READ");
505 MUONLoader -> LoadRecPoints("READ");
508 //--------------------------------------------------------------------------------------------------------------
509 //Set mag field; waiting for mag field in CDB
510 printf("Loading field map...\n");
511 if (!AliTracker::GetFieldMap()) {
512 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
513 AliTracker::SetFieldMap(field, kFALSE);}
516 //--------------------------------------------------------------------------------------------------------------
517 //Set Field Map for track extrapolation
518 AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
521 //Creating a MUON data container
522 //--------------------------------------------------------------------------------------------------------------
523 AliMUONData muondata(MUONLoader,"MUON","MUON");
526 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
529 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
531 if ( event2Check!=0 )
533 printf("\r>>> Event %d ",iEvent);
534 RunLoader->GetEvent(iEvent);
537 muondata.SetTreeAddress("RT");
538 muondata.GetRecTracks();
539 tracks = muondata.RecTracks();
542 nTracks = (Int_t) tracks -> GetEntriesFast();
543 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
545 track = (AliMUONTrack*) tracks ->At(iTrack) ;
546 trackParams = track ->GetTrackParamAtHit();
547 hits = track ->GetHitForRecAtHit() ;
548 chamber = firstChamber - 1;
553 nHits = (Int_t) hits -> GetEntriesFast();
554 for ( iHit = 0; iHit < nHits; ++iHit )
556 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
557 hit = (AliMUONHitForRec*) hits -> At(iHit);
558 chamberNbr = hit -> GetChamberNumber();
560 if ( chamberNbr == oldChamber )
563 oldChamber = chamberNbr;
564 if ( chamberNbr > chamber - k ) {
565 if ( chamber < lastChamber ) {
566 if ( chamberNbr == chamber ) {
568 Double_t Px, Py, Pz, Pr;
574 Px = trackParam -> Px() ;
575 Py = trackParam -> Py() ;
576 Pz = trackParam -> Pz() ;
577 Pr = sqrt( Px*Px + Py*Py );
579 //The calculation of theta, the angle of incidence:
580 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
582 if ( theta < 79 ) iTheta = 11;
583 else if ( theta < 90 ) iTheta = int( theta - 79.);
586 eff [iTheta] = eff [iTheta] + 1;
587 trackN [iTheta] = trackN [iTheta] + 1;
588 chamber = chamber + 1;
593 Double_t chamberZpos;
596 Double_t Px, Py, Pz, Pr;
602 chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
603 AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
605 Px = trackParam -> Px() ;
606 Py = trackParam -> Py() ;
607 Pz = trackParam -> Pz() ;
608 Pr = sqrt( Px*Px + Py*Py );
610 //The calculation of thetaI, the angle of incidence:
611 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
613 if ( theta < 79 ) iTheta = 11;
614 else if ( theta < 90 ) iTheta = int( theta - 79.);
617 eff [iTheta] = eff [iTheta] + 0;
618 trackN [iTheta] = trackN [iTheta] + 1;
619 chamber = chamber + 1;
626 if ( iHit == nHits-1 ) {
627 if ( chamber == lastChamber )
641 muondata.ResetRecTracks();
642 if ( event2Check!=0 )
648 Double_t efficiency [11];
652 TGraph * effTheta = new TGraph ();
654 if ( firstChamber == lastChamber ) {
655 for ( th = 0; th < 11; ++th )
657 cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
658 cout<<"The chamber "<<firstChamber<<" has responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
661 Double_t e = eff [th];
662 Double_t n = trackN [th];
665 efficiency [th] = 0.;
666 cout<<"Efficiency = ? % (IS UNKNOWN) \n";
669 efficiency [th] = 100.*e/n;
670 cout<<"Efficiency = "<<efficiency [th]<<" %\n";
676 for ( th = 0; th < 11; ++th )
678 cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
679 cout<<"The chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
682 Double_t e = eff [th];
683 Double_t n = trackN [th];
686 efficiency [th] = 0.;
687 cout<<"Efficiency = ? % (IS UNKNOWN) \n";
690 efficiency [th] = 100.*e/n;
691 cout<<"Efficiency = "<<efficiency [th]<<" %\n";
696 TCanvas * c1 = new TCanvas ();
697 effTheta = new TGraph( i, t, efficiency );
699 effTheta -> Draw("ALP");
701 MUONLoader->UnloadTracks();
710 /*****************************************************************************************************************/
711 /*****************************************************************************************************************/
714 void resolution( Int_t event2Check=0, char * filename="galice.root" )
716 cout<<"\nChamber(s) chosen;\nFirst chamber:";
718 cout<<"Last chamber:";
726 hDeltax = new TH1F ("hDeltax " , "No Interest", 100, -0.3, 0.3);
727 hDeltay = new TH1F ("hDeltay " , "No Interest", 500, -0.1, 0.1);
728 hDelta3D = new TH2F ("hDelta3D" , "No Interest", 100, -0.3, 0.3, 500, -0.1, 0.1);
731 //Creating Run Loader and openning file containing Hits
732 //--------------------------------------------------------------------------------------------------------------
733 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
734 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
738 //--------------------------------------------------------------------------------------------------------------
739 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
740 MUONLoader -> LoadTracks("READ");
741 MUONLoader -> LoadRecPoints("READ");
744 //Creating a MUON data container
745 //--------------------------------------------------------------------------------------------------------------
746 AliMUONData muondata(MUONLoader,"MUON","MUON");
749 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
752 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
756 printf("\r>>> Event %d ",iEvent);
757 RunLoader->GetEvent(iEvent);
760 muondata.SetTreeAddress("RT");
761 muondata.GetRecTracks();
762 tracks = muondata.RecTracks();
765 nTracks = (Int_t) tracks -> GetEntriesFast();
766 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
768 track = (AliMUONTrack*) tracks ->At(iTrack) ;
769 trackParams = track ->GetTrackParamAtHit();
770 hits = track ->GetHitForRecAtHit() ;
773 nHits = (Int_t) hits -> GetEntriesFast();
774 for ( iHit = 0; iHit < nHits; ++iHit )
776 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
777 hit = (AliMUONHitForRec*) hits -> At(iHit);
778 chamberNbr = hit -> GetChamberNumber();
779 if ( chamberNbr >= firstChamber-1 ) {
780 if ( chamberNbr < lastChamber ) {
786 Double_t deltaX, deltaY;
788 traX = trackParam -> GetNonBendingCoor();
789 traY = trackParam -> GetBendingCoor();
790 hitX = hit -> GetNonBendingCoor();
791 hitY = hit -> GetBendingCoor();
793 deltaX = traX - hitX;
794 deltaY = traY - hitY;
796 hDeltax -> Fill (deltaX);
797 hDeltay -> Fill (deltaY);
798 hDelta3D -> Fill (deltaX, deltaY);
805 muondata.ResetRecTracks();
806 if ( event2Check!=0 )
815 if ( firstChamber == lastChamber ) {
816 sprintf(hXtitle,"DeltaX Chamber %d",firstChamber);
817 sprintf(hYtitle,"DeltaY Chamber %d",lastChamber);
818 sprintf(h3title,"DeltaX and Delta Y Chamber %d",lastChamber);
821 sprintf(hXtitle,"DeltaX Chamber %d to %d",firstChamber,lastChamber);
822 sprintf(hYtitle,"DeltaY Chamber %d to %d",firstChamber,lastChamber);
823 sprintf(h3title,"DeltaX and Delta Y Chamber %d to %d",firstChamber,lastChamber);
827 hDeltax -> SetTitle (hXtitle);
828 hDeltay -> SetTitle (hYtitle);
829 hDelta3D -> SetTitle (h3title);
831 Double_t rmsX = hDeltax -> GetRMS ();
832 Double_t errX = hDeltax -> GetRMSError();
834 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
835 hDeltay -> Fit("fY","R","E");
836 Double_t sigY = fY -> GetParameter(2);
837 Double_t errY = fY -> GetParError (2);
839 if ( firstChamber == lastChamber ) {
840 cout<<"\nThe resolution (Non Bending plane) of the chamber "<<firstChamber<<" = "<<rmsX<<" cm +-"<<errX;
841 cout<<"\nThe resolution (Bending plane) of the chamber "<<firstChamber<<" = "<<sigY<<" cm +- "<<errY;
844 cout<<"\nThe resolution (Non Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<rmsX<<" cm +- "<<errX;
845 cout<<"\nThe resolution (Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<sigY<<" cm +- "<<errY;
849 MUONLoader->UnloadTracks();
859 /*****************************************************************************************************************/
860 /*****************************************************************************************************************/
861 /*RESOLUTION IN FUNCTION OF PHI*/
863 void resolutionPhi( Int_t event2Check=0, char * filename="galice.root" )
865 cout<<"\nChamber(s) chosen;\nFirst chamber:";
867 cout<<"Last chamber:";
874 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
875 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
876 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
877 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
878 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
880 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
881 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
882 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
883 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
884 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
887 //Creating Run Loader and openning file containing Hits
888 //--------------------------------------------------------------------------------------------------------------
889 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
890 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
894 //--------------------------------------------------------------------------------------------------------------
895 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
896 MUONLoader -> LoadTracks("READ");
897 MUONLoader -> LoadRecPoints("READ");
900 //Creating a MUON data container
901 //--------------------------------------------------------------------------------------------------------------
902 AliMUONData muondata(MUONLoader,"MUON","MUON");
905 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
908 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
910 if ( event2Check!=0 )
912 printf("\r>>> Event %d ",iEvent);
913 RunLoader->GetEvent(iEvent);
916 muondata.SetTreeAddress("RT");
917 muondata.GetRecTracks();
918 tracks = muondata.RecTracks();
921 nTracks = (Int_t) tracks -> GetEntriesFast();
922 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
924 track = (AliMUONTrack*) tracks ->At(iTrack) ;
925 trackParams = track ->GetTrackParamAtHit();
926 hits = track ->GetHitForRecAtHit() ;
929 nHits = (Int_t) hits -> GetEntriesFast();
930 for ( iHit = 0; iHit < nHits; ++iHit )
932 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
933 hit = (AliMUONHitForRec* ) hits -> At(iHit);
934 chamberNbr = hit -> GetChamberNumber();
935 if ( chamberNbr >= firstChamber -1 ) {
936 if ( chamberNbr < lastChamber ) {
950 traX = trackParam -> GetNonBendingCoor();
951 traY = trackParam -> GetBendingCoor() ;
953 hitX = hit -> GetNonBendingCoor();
954 hitY = hit -> GetBendingCoor() ;
956 aveX = 1./2.*(traX + hitX);
957 aveY = 1./2.*(traY + hitY);
959 //The calculation of phi:
960 phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));
962 if ( phi < 0.) phi = 360. - fabs(phi);
963 iPhi = int( phi/72. );
965 deltaX = traX - hitX;
966 deltaY = traY - hitY;
968 hDeltaX [iPhi] -> Fill( deltaX );
969 hDeltaY [iPhi] -> Fill( deltaY );
978 muondata.ResetRecTracks();
979 if ( event2Check!=0 )
988 Int_t phiMin, phiMax;
997 for ( iPhi = 0; iPhi < iPhiMax; ++iPhi )
1003 phiMax = 72*iPhi + 72;
1005 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1007 if ( firstChamber == lastChamber ) {
1008 sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d",phiMin,phiMax,firstChamber);
1009 sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d",phiMin,phiMax,lastChamber);
1012 sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
1013 sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
1016 hDeltaY [iPhi] -> SetTitle(hYtitle);
1017 hDeltaX [iPhi] -> SetTitle(hXtitle);
1019 hDeltaY [iPhi] -> Fit("fY","R","E");
1020 sigmaY [iPhi] = fY -> GetParameter(2) ;
1021 errSY [iPhi] = fY -> GetParError(2) ;
1023 rmsX [iPhi] = hDeltaX [iPhi] -> GetRMS();
1024 errSX [iPhi] = hDeltaX [iPhi] -> GetRMSError();
1026 phi [iPhi] = 72*iPhi + 36 ;
1030 //--------------------------------------------------------------------------------------------------------------
1031 //For plotting resolution in function of the angle of incidence
1033 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1037 TGraphErrors * GraphX = new TGraphErrors( iPhiMax, phi, rmsX, errPh, errSX);
1038 GraphX->SetTitle("Resolution en X (Phi)");
1039 GraphX->Draw("ALP");
1042 TGraphErrors * GraphY = new TGraphErrors( iPhiMax, phi, sigmaY, errPh, errSY);
1043 GraphY->SetTitle("Resolution en Y (Phi)");
1044 GraphY->Draw("ALP");
1048 MUONLoader->UnloadTracks();
1058 /*****************************************************************************************************************/
1059 /*****************************************************************************************************************/
1060 /*RESOLUTION IN FUNCTION OF THETA*/
1062 void resolutionTheta( Int_t event2Check=0, char * filename="galice.root" )
1064 cout<<"\nChamber(s) chosen;\nFirst chamber:";
1066 cout<<"Last chamber:";
1073 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
1074 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
1075 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
1076 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
1077 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
1078 hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
1079 hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
1080 hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
1081 hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
1083 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
1084 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
1085 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
1086 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
1087 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
1088 hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
1089 hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
1090 hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
1091 hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
1094 //Creating Run Loader and openning file containing Hits
1095 //--------------------------------------------------------------------------------------------------------------
1096 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
1097 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
1100 //Getting MUONLoader
1101 //--------------------------------------------------------------------------------------------------------------
1102 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1103 MUONLoader -> LoadTracks("READ");
1104 MUONLoader -> LoadRecPoints("READ");
1107 //Creating a MUON data container
1108 //--------------------------------------------------------------------------------------------------------------
1109 AliMUONData muondata(MUONLoader,"MUON","MUON");
1112 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1115 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1117 if ( event2Check!=0 )
1119 printf("\r>>> Event %d ",iEvent);
1120 RunLoader->GetEvent(iEvent);
1123 muondata.SetTreeAddress("RT");
1124 muondata.GetRecTracks();
1125 tracks = muondata.RecTracks();
1128 nTracks = (Int_t) tracks -> GetEntriesFast();
1129 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1131 track = (AliMUONTrack*) tracks ->At(iTrack) ;
1132 trackParams = track ->GetTrackParamAtHit();
1133 hits = track ->GetHitForRecAtHit() ;
1136 nHits = (Int_t) hits -> GetEntriesFast();
1137 for ( iHit = 0; iHit < nHits; ++iHit )
1139 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
1140 hit = (AliMUONHitForRec* ) hits -> At(iHit);
1141 chamberNbr = hit -> GetChamberNumber();
1142 if ( chamberNbr >= firstChamber -1 ) {
1143 if ( chamberNbr < lastChamber ) {
1145 Double_t traX, traY;
1146 Double_t hitX, hitY, hitZ;
1147 Double_t aveY, aveX;
1157 traX = trackParam -> GetNonBendingCoor();
1158 traY = trackParam -> GetBendingCoor() ;
1160 hitX = hit -> GetNonBendingCoor();
1161 hitY = hit -> GetBendingCoor() ;
1162 hitZ = hit -> GetZ();
1164 aveX = 1./2.*(traX + hitX);
1165 aveY = 1./2.*(traY + hitY);
1167 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
1168 theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ ));
1170 iTheta = int( theta );
1171 if( theta > 10 ) iTheta = 9;
1172 if( theta < 1 ) iTheta = 1;
1174 deltaX = traX - hitX;
1175 deltaY = traY - hitY;
1177 hDeltaX [iTheta-1] -> Fill( deltaX );
1178 hDeltaY [iTheta-1] -> Fill( deltaY );
1186 //End loop on tracks
1187 muondata.ResetRecTracks();
1188 if ( event2Check!=0 )
1192 //End loop on events
1196 Int_t iThetaMax = 9;
1197 Int_t thetaMin, thetaMax;
1206 for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1211 //To find the polar angle
1212 thetaMin = 178 - iTheta;
1213 thetaMax = 179 - iTheta;
1215 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1217 if ( firstChamber == lastChamber ) {
1218 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
1219 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
1222 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1223 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1226 hDeltaY [iTheta] -> SetTitle(hYtitle);
1227 hDeltaX [iTheta] -> SetTitle(hXtitle);
1229 hDeltaY [iTheta] -> Fit("fY","R","E");
1230 sigmaY [iTheta] = fY -> GetParameter(2) ;
1231 errSY [iTheta] = fY -> GetParError(2) ;
1233 rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS();
1234 errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1236 theta [iTheta] = 178.5 - iTheta;
1237 ErrTh [iTheta] = 0.5;
1240 //--------------------------------------------------------------------------------------------------------------
1241 //For plotting resolution in function of the angle of incidence
1243 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1247 TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, ErrTh, errSX);
1248 GraphX->SetTitle("Resolution en X (Theta)");
1249 GraphX->Draw("ALP");
1252 TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, ErrTh, errSY);
1253 GraphY->SetTitle("Resolution en Y (Theta)");
1254 GraphY->Draw("ALP");
1257 MUONLoader->UnloadTracks();
1265 /*****************************************************************************************************************/
1266 /*****************************************************************************************************************/
1267 /*RESOLUTION IN FUNCTION OF THETA OF INCIDENCE*/
1269 void resolutionThetaI( Int_t event2Check=0, char * filename="galice.root" )
1271 cout<<"\nChamber(s) chosen;\nFirst chamber:";
1273 cout<<"Last chamber:";
1280 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
1281 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
1282 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
1283 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
1284 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
1285 hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
1286 hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
1287 hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
1288 hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
1289 hDeltaX[9] = new TH1F ("hDeltaX9" , "No Interest", 100, -0.3, 0.3);
1290 hDeltaX[10]= new TH1F ("hDeltaX10", "No Interest", 100, -0.3, 0.3);
1292 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
1293 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
1294 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
1295 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
1296 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
1297 hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
1298 hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
1299 hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
1300 hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
1301 hDeltaY[9] = new TH1F ("hDeltaY9" , "No Interest", 500, -0.1, 0.1);
1302 hDeltaY[10]= new TH1F ("hDeltaY10", "No Interest", 500, -0.1, 0.1);
1305 //Creating Run Loader and openning file containing Hits
1306 //--------------------------------------------------------------------------------------------------------------
1307 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
1308 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
1311 //Getting MUONLoader
1312 //--------------------------------------------------------------------------------------------------------------
1313 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1314 MUONLoader -> LoadTracks("READ");
1315 MUONLoader -> LoadRecPoints("READ");
1318 //Creating a MUON data container
1319 //--------------------------------------------------------------------------------------------------------------
1320 AliMUONData muondata(MUONLoader,"MUON","MUON");
1323 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1326 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1328 if ( event2Check!=0 )
1330 printf("\r>>> Event %d ",iEvent);
1331 RunLoader->GetEvent(iEvent);
1334 muondata.SetTreeAddress("RT");
1335 muondata.GetRecTracks();
1336 tracks = muondata.RecTracks();
1339 nTracks = (Int_t) tracks -> GetEntriesFast();
1340 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1342 track = (AliMUONTrack*) tracks ->At(iTrack) ;
1343 trackParams = track ->GetTrackParamAtHit();
1344 hits = track ->GetHitForRecAtHit() ;
1347 nHits = (Int_t) hits -> GetEntriesFast();
1348 for ( iHit = 0; iHit < nHits; ++iHit )
1350 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
1351 hit = (AliMUONHitForRec* ) hits -> At(iHit);
1352 chamberNbr = hit -> GetChamberNumber();
1353 if ( chamberNbr >= firstChamber - 1 ) {
1354 if ( chamberNbr < lastChamber ) {
1356 Double_t Px, Py, Pz, Pr;
1359 Double_t traX, traY;
1360 Double_t hitX, hitY;
1370 Px = trackParam -> Px() ;
1371 Py = trackParam -> Py() ;
1372 Pz = trackParam -> Pz() ;
1373 Pr = sqrt( Px*Px + Py*Py );
1375 traX = trackParam -> GetNonBendingCoor();
1376 traY = trackParam -> GetBendingCoor();
1378 hitX = hit -> GetNonBendingCoor();
1379 hitY = hit -> GetBendingCoor();
1381 //The calculation of theta, the angle of incidence:
1382 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr))
1384 if ( theta < 79 ) iTheta = 0;
1385 else iTheta = int( theta - 79.);
1387 deltaX = traX - hitX;
1388 deltaY = traY - hitY;
1390 hDeltaX [iTheta] -> Fill (deltaX);
1391 hDeltaY [iTheta] -> Fill (deltaY);
1399 //End loop on tracks
1400 muondata.ResetRecTracks();
1401 if ( event2Check!=0 )
1405 //End loop on events
1409 Int_t iThetaMax = 11;
1410 Int_t thetaMin, thetaMax;
1419 for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1424 thetaMin = iTheta + 79;
1425 thetaMax = iTheta + 80;
1427 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1429 if ( firstChamber == lastChamber ) {
1430 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
1431 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
1434 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1435 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1438 hDeltaY [iTheta] -> SetTitle(hYtitle);
1439 hDeltaX [iTheta] -> SetTitle(hXtitle);
1441 hDeltaY [iTheta] -> Fit("fY","R","E");
1442 sigmaY [iTheta] = fY -> GetParameter(2) ;
1443 errSY [iTheta] = fY -> GetParError(2) ;
1445 rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS();
1446 errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1448 theta [iTheta] = iTheta + 79.5 ;
1449 errTh [iTheta] = 0.5;
1452 //--------------------------------------------------------------------------------------------------------------
1453 //For plotting resolution in function of the angle of incidence
1455 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1459 TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, errTh, errSX);
1460 GraphX->SetTitle("Resolution en X (Theta)");
1461 GraphX->Draw("ALP");
1464 TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, errTh, errSY);
1465 GraphY->SetTitle("Resolution en Y (Theta)");
1466 GraphY->Draw("ALP");
1469 MUONLoader->UnloadTracks();