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