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