Macro to calculate the resolution and the efficiency of chamber(s) (Nicolas)
authormartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 29 Jan 2007 10:03:58 +0000 (10:03 +0000)
committermartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 29 Jan 2007 10:03:58 +0000 (10:03 +0000)
MUON/MUONResoEffChamber.C [new file with mode: 0644]

diff --git a/MUON/MUONResoEffChamber.C b/MUON/MUONResoEffChamber.C
new file mode 100644 (file)
index 0000000..4234e73
--- /dev/null
@@ -0,0 +1,1458 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+// ROOT includes
+#include "TBranch.h"
+#include "TClonesArray.h"
+#include "TTree.h"
+#include "TNtuple.h"
+#include "TParticle.h"
+#include "TFile.h"
+
+#include "TH1.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TMath.h"
+
+#include "TCanvas.h"
+#include "TGraph.h"
+#include "TGraphErrors.h"
+#include "TGraph2D.h"
+#include "TStyle.h"
+#include "TFitter.h"
+#include "TRandom.h"
+
+// STEER includes
+#include "AliRun.h"
+#include "AliRunLoader.h"
+#include "AliHeader.h"
+#include "AliLoader.h"
+#include "AliTracker.h"
+#include "AliStack.h"
+#include "AliMagFMaps.h"
+
+
+// MUON includes
+#include "AliMUON.h"
+#include "AliMUONData.h"
+#include "AliMUONConstants.h"
+
+#include "AliMUONHit.h"
+#include "AliMUONHitForRec.h"
+#include "AliMUONRawCluster.h"
+#include "AliMUONTrack.h"
+#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
+
+#include "AliMpVSegmentation.h"
+#include "AliMpIntPair.h"
+#include "AliMpDEManager.h"
+#endif
+
+
+//Macro to calculate the resolution and the efficiency of chamber(s) chosen in function
+//of Phi (angle on the chamber between the X axis and the straight line created by the
+//center of the chamber and the impact of particles on this chamber, the azimuthal angle)
+//and Theta (the polar angle), or in function of ThetaI (angle of incidence of particles
+//on the chamber)
+
+
+
+
+  const Double_t kInvPi = 1./3.14159265;
+
+
+//Chamber number:
+  Int_t chamberNbr;
+//Number of events:
+  Int_t nEvents, iEvent;
+//Number of tracks:
+  Int_t nTracks, iTrack;
+//Number of hits:
+  Int_t nHits,iHit;
+//Chamber(s) chosen
+  Int_t firstChamber, lastChamber;
+
+  AliMUONTrack      * track  ;
+  AliMUONHitForRec  * hit = 0;
+  AliMUONTrackParam * trackParam = 0;
+  TClonesArray      * tracks  ;
+  TClonesArray      * trackParams;
+  TClonesArray      * hits  ;
+
+
+
+
+
+
+/*****************************************************************************************************************/
+/*****************************************************************************************************************/
+                                               /*EFFICIENCY*/
+
+void efficiency( Int_t event2Check=0, char * filename="galice.root" )
+{
+  cout<<"\nChamber(s) chosen;\nFirst chamber:";
+  cin>>firstChamber;
+  cout<<"Last chamber:";
+  cin>>lastChamber;
+  cout<<"\n\n";
+
+//Creating Run Loader and openning file containing Hits
+//--------------------------------------------------------------------------------------------------------------
+  AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
+  if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
+
+
+//Getting MUONLoader
+//--------------------------------------------------------------------------------------------------------------
+  AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
+  MUONLoader -> LoadTracks("READ");
+  MUONLoader -> LoadRecPoints("READ");
+
+
+//Creating a MUON data container
+//--------------------------------------------------------------------------------------------------------------
+  AliMUONData muondata(MUONLoader,"MUON","MUON");
+
+
+  nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
+
+    //Loop on events
+    Int_t trackNb = 0;
+    Int_t chamber[10] = {0};
+    Int_t detEltNew, detElt;
+    Int_t detEltOld = 0;
+
+    for ( iEvent = 0; iEvent < nEvents; ++iEvent )
+    {
+      if ( event2Check!=0 ) 
+         iEvent=event2Check;
+      printf("\r>>> Event %d ",iEvent);
+      RunLoader -> GetEvent(iEvent);
+      //Addressing
+      muondata.SetTreeAddress("RT");
+      muondata.GetRecTracks();
+      tracks = muondata.RecTracks();
+
+      //Loop on track
+      nTracks = (Int_t) tracks -> GetEntriesFast();
+      for ( iTrack = 0; iTrack < nTracks; ++iTrack )
+      { 
+       track     = (AliMUONTrack*) tracks -> At(iTrack);
+       hits = track -> GetHitForRecAtHit(); 
+       detEltOld = 0; 
+       //Loop on hit
+       nHits = (Int_t) hits -> GetEntriesFast();
+
+       for ( iHit = 0; iHit < nHits; ++iHit )
+       { 
+         hit        = (AliMUONHitForRec*) hits -> At(iHit);
+         chamberNbr =                     hit  -> GetChamberNumber();
+         detElt     =                     hit  -> GetDetElemId();
+         detEltNew  = int(detElt/100);
+         if( chamberNbr >= firstChamber-1 ) {
+           if( chamberNbr < lastChamber ) {
+             if( detEltNew == detEltOld ) 
+                 continue;
+             else {
+               chamber[chamberNbr] = chamber[chamberNbr] + 1; 
+               detEltOld = detEltNew;
+             }
+           }
+         }
+       }           
+       //End loop on hit
+
+      }
+      //End loop on track 
+      muondata.ResetRecTracks();
+      if (event2Check != 0)
+         iEvent = nEvents;
+      trackNb = trackNb + nTracks;
+    }
+    //End loop on event
+//--------------------------------------------------------------------------------------------------------------
+
+    cout<<"\n\n\n";
+    for (Int_t i = firstChamber-1; i < lastChamber; ++i )
+    { 
+      printf ("\nChamber %d has responded %d times on %d tracks", i+1, chamber[i], trackNb);
+      if (trackNb == 0) 
+         cout<<"\nEfficiency = ? (IS UNKNOWN)\n";
+      else {
+       Double_t eff = chamber[i]*100./trackNb; cout<<"\nEfficiency = "<<eff<<" %\n";
+      }
+    }
+    cout<<"\n\n\n";
+    MUONLoader->UnloadTracks();
+}
+
+
+
+
+
+/*****************************************************************************************************************/
+/*****************************************************************************************************************/
+                               /*EFFICIENCY IN FUNCTION OF THETA AND PHI*/
+
+void efficiencyThetaPhi( Int_t event2Check=0, char * filename="galice.root" )
+{
+  cout<<"\nChamber(s) chosen;\nFirst chamber:";
+  cin>>firstChamber;
+  cout<<"Last chamber:";
+  cin>>lastChamber;
+  cout<<"\n\n";
+
+  Int_t eff    [54] = {0};
+  Int_t trackN [54] = {0};
+  Int_t chamber;
+  Int_t oldChamber;
+
+//Creating Run Loader and openning file containing Hits
+//--------------------------------------------------------------------------------------------------------------
+  AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
+  if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
+
+
+//Getting MUONLoader
+//--------------------------------------------------------------------------------------------------------------
+  AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
+  MUONLoader -> LoadTracks("READ");
+  MUONLoader -> LoadRecPoints("READ");
+
+
+//--------------------------------------------------------------------------------------------------------------
+//Set mag field; waiting for mag field in CDB 
+  printf("Loading field map...\n");
+  if (!AliTracker::GetFieldMap()) {
+  AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
+  AliTracker::SetFieldMap(field, kFALSE);}
+
+
+//--------------------------------------------------------------------------------------------------------------
+//Set Field Map for track extrapolation
+  AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
+
+
+//Creating a MUON data container
+//--------------------------------------------------------------------------------------------------------------
+  AliMUONData muondata(MUONLoader,"MUON","MUON");
+
+
+  nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
+
+    //Loop on events
+    for ( iEvent = 0; iEvent < nEvents; ++iEvent )
+    { 
+      if ( event2Check!=0 )
+         iEvent=event2Check;
+      printf("\r>>> Event %d ",iEvent);
+      RunLoader->GetEvent(iEvent);
+      
+      //Addressing
+      muondata.SetTreeAddress("RT");     
+      muondata.GetRecTracks();
+      tracks = muondata.RecTracks();
+
+      //Loop on track
+      nTracks = (Int_t) tracks -> GetEntriesFast();
+      for ( iTrack = 0; iTrack < nTracks; ++iTrack )
+      {
+       track       = (AliMUONTrack*) tracks ->At(iTrack)          ;
+       trackParams =                 track  ->GetTrackParamAtHit();
+       hits        =                 track  ->GetHitForRecAtHit() ;
+       chamber     = firstChamber-1;
+       oldChamber  = -1;
+       Int_t k     =  1;
+  
+       //Loop on hits
+       nHits = (Int_t) hits->GetEntriesFast();
+       for ( iHit = 0; iHit<nHits; ++iHit )
+       { 
+         trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
+         hit        = (AliMUONHitForRec* ) hits        -> At(iHit);
+         chamberNbr = hit -> GetChamberNumber();
+              
+         if ( chamberNbr == oldChamber ) 
+             continue;
+         else { 
+           oldChamber = chamberNbr;
+           if ( chamberNbr > chamber - k ) {
+             if ( chamber < lastChamber ) {
+               if ( chamberNbr == chamber ) {
+                 //Positions
+                 Double_t traX, traY, traZ;
+                 Double_t hitX, hitY, hitZ;
+                 Double_t aveX, aveY      ;
+
+                 //Angle (Phi)
+                 Double_t phi   = 0.;
+                 Double_t theta = 0.;
+                 Int_t    iPhi   = 0 ;
+                 Int_t    iTheta = 0 ;                                      
+                                                 
+                 traX = trackParam -> GetNonBendingCoor();
+                 traY = trackParam -> GetBendingCoor()   ;
+                 traZ = trackParam -> GetZ()             ;
+                                                 
+                 hitX = hit        -> GetNonBendingCoor();
+                 hitY = hit        -> GetBendingCoor()   ;
+                 hitZ = hit        -> GetZ()             ;
+                                                 
+                 aveX = 1./2.*(traX + hitX);
+                 aveY = 1./2.*(traY + hitY);
+               
+                 //The calculation of phi:
+                 phi   = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));
+
+                 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
+                 theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ ));
+
+                 if ( phi < 0.) phi = 360 - abs(phi);
+                 iPhi = int( phi/72. );
+                 iTheta = int( theta );
+                 if( theta > 10 ) iTheta = 9;
+                 if( theta < 1  ) iTheta = 1;
+                                              
+                 eff    [9*iPhi+iTheta-1] = eff    [9*iPhi+iTheta-1] + 1;
+                 trackN [9*iPhi+iTheta-1] = trackN [9*iPhi+iTheta-1] + 1;  
+                 chamber                  = chamber + 1;
+               }   
+
+               else {
+                 //Positions
+                 Double_t chamberZpos;
+                 Double_t exXpos, exYpos;
+                                                 
+                 //Angles
+                 Double_t phi   = 0.;
+                 Double_t theta = 0.;
+                 Int_t    iPhi   = 0 ;
+                 Int_t    iTheta = 0 ;
+                                                 
+                 chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
+                 AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
+                 exXpos      = (Double_t) trackParam->GetNonBendingCoor();
+                 exYpos      = (Double_t) trackParam->GetBendingCoor();
+                 
+                 //The calculation of phi:
+                 phi   = 180. * kInvPi * (TMath::ATan2( exYpos, exXpos ));                               
+
+                 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
+                 theta = 180. * kInvPi * (TMath::ATan2( sqrt(exXpos*exXpos+exYpos*exYpos), - chamberZpos ));
+
+                 if ( phi < 0.) phi = 360 - abs(phi);
+                 iPhi = int( phi/72. );
+                 iTheta = int( theta );
+                 if( theta > 10 ) iTheta = 9;
+                 if( theta < 1  ) iTheta = 1;
+
+                 eff    [9*iPhi+iTheta-1] = eff    [9*iPhi+iTheta-1] + 0;
+                 trackN [9*iPhi+iTheta-1] = trackN [9*iPhi+iTheta-1] + 1;
+                 chamber                  = chamber + 1;
+                 iHit                     = iHit - 1;
+               }
+             }
+           }
+         }
+                                                  
+         if ( iHit == nHits-1 ) {
+           if    ( chamber == lastChamber )
+               continue;
+           else {
+             oldChamber = -1;
+             k          =  5;
+             iHit       = iHit-1;
+           }
+         }                                   
+       } 
+       //End Loop on hits
+
+      } 
+      //End Loop on tracks
+
+      muondata.ResetRecTracks();
+      if ( event2Check!=0 )
+         iEvent=nEvents;
+    }
+    //End Loop on events
+
+
+    TCanvas * c1 = new TCanvas();
+    TGraph2D* effPhiTheta = new TGraph2D();
+    Double_t  efficiency = 0;
+
+    if ( firstChamber == lastChamber )
+    {
+      for ( Int_t ph = 0; ph < 5; ++ph )
+      {
+       for ( Int_t th = 1; th < 10; ++th )
+       {
+         Int_t i = 9*ph+th-1;
+         cout<<"\nFor Phi = ["<<ph*72<<","<<ph*72+72<<"] and Theta = ["<<179-th<<","<<180-th<<"]:";
+         cout<<"\nThe chamber "<<firstChamber<<" has responded "<<eff[i]<<" times on "<<trackN[i]<<" tracks\n";
+
+         Double_t e = eff   [i] ;
+         Double_t n = trackN[i] ;
+         Double_t p = ph*72.+36.;
+         Double_t t = th*1. +0.5;
+       
+         if ( trackN[i] == 0 ) {
+           efficiency = 0.;
+           cout<<"Efficiency = ? % ( IS UNKNOWN )\n";
+         }
+         else {
+           efficiency = 100.*e/n;
+           cout<<"Efficiency = "<<efficiency<<" %\n";
+         }
+
+         effPhiTheta -> SetPoint( i, p, t, efficiency);
+       }
+      }
+    }
+
+    else{
+      for ( Int_t ph = 0; ph < 5; ++ph )
+      {
+       for ( Int_t th = 1; th < 10; ++th )
+       {
+         Int_t i = 9*ph+th-1;
+         cout<<"\nFor Phi = ["<<ph*72<<","<<ph*72+72<<"] and Theta = ["<<179-th<<","<<180-th<<"]:";
+         cout<<"\nThe chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff[i]<<" times on "<<trackN[i]<<" tracks\n";
+
+         Double_t e = eff   [i] ;
+         Double_t n = trackN[i] ;
+         Double_t p = ph*72.+36.;
+         Double_t t = th*1. +0.5;
+       
+         if ( trackN[i] == 0 ) {
+           efficiency = 0.;
+           cout<<"Efficiency = ? % ( IS UNKNOWN )\n";
+         }
+         else {
+           efficiency = 100.*e/n;
+           cout<<"Efficiency = "<<efficiency<<" %\n";
+         }
+         
+         effPhiTheta -> SetPoint( i, p, t, efficiency);
+       }
+      }
+    }
+
+    gStyle->SetPalette(1);
+    effPhiTheta -> Draw("surf1");
+
+    cout<<"\n\n\n";
+    MUONLoader->UnloadTracks();
+    c1->Update();
+
+}
+
+
+
+
+
+
+
+/*****************************************************************************************************************/
+/*****************************************************************************************************************/
+                          /*EFFICIENCY IN FUNCTION OF THETA OF INCIDENCE*/
+
+
+void efficiencyThetaI( Int_t event2Check=0, char * filename="galice.root" )
+{
+  cout<<"\nChamber(s) chosen;\nFirst chamber:";
+  cin>>firstChamber;
+  cout<<"Last chamber:";
+  cin>>lastChamber;
+  cout<<"\n\n";
+
+  Int_t eff    [12] = {0};
+  Int_t trackN [12] = {0};
+  Int_t chamber;
+  Int_t oldChamber;
+
+
+//Creating Run Loader and openning file containing Hits
+//--------------------------------------------------------------------------------------------------------------
+  AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
+  if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
+
+
+//Getting MUONLoader
+//--------------------------------------------------------------------------------------------------------------
+  AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
+  MUONLoader -> LoadTracks("READ");
+  MUONLoader -> LoadRecPoints("READ");
+
+
+//--------------------------------------------------------------------------------------------------------------
+//Set mag field; waiting for mag field in CDB 
+  printf("Loading field map...\n");
+  if (!AliTracker::GetFieldMap()) {
+  AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
+  AliTracker::SetFieldMap(field, kFALSE);}
+
+
+//--------------------------------------------------------------------------------------------------------------
+//Set Field Map for track extrapolation
+  AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
+
+
+//Creating a MUON data container
+//--------------------------------------------------------------------------------------------------------------
+  AliMUONData muondata(MUONLoader,"MUON","MUON");
+
+
+  nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
+
+    //Loop on events
+    for ( iEvent = 0; iEvent < nEvents; ++iEvent )
+    { 
+      if ( event2Check!=0 )
+         iEvent=event2Check;
+      printf("\r>>> Event %d ",iEvent);
+      RunLoader->GetEvent(iEvent);
+      
+      //Addressing
+      muondata.SetTreeAddress("RT");     
+      muondata.GetRecTracks();
+      tracks = muondata.RecTracks();
+
+      //Loop on track
+      nTracks = (Int_t) tracks -> GetEntriesFast();
+      for ( iTrack = 0; iTrack < nTracks; ++iTrack )
+      { 
+       track       = (AliMUONTrack*) tracks ->At(iTrack)          ;
+       trackParams =                 track  ->GetTrackParamAtHit();
+       hits        =                 track  ->GetHitForRecAtHit() ;
+       chamber     = firstChamber - 1;
+       oldChamber  = -1;
+       Int_t k     =  1;
+
+       //Loop on hits
+       nHits = (Int_t) hits -> GetEntriesFast();
+       for ( iHit = 0; iHit < nHits; ++iHit )
+       { 
+         trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
+         hit        = (AliMUONHitForRec*)  hits        -> At(iHit);
+         chamberNbr = hit -> GetChamberNumber();
+
+         if ( chamberNbr == oldChamber )
+             continue;
+         else {
+           oldChamber = chamberNbr;
+           if ( chamberNbr > chamber - k ) {
+             if ( chamber < lastChamber ) {
+               if ( chamberNbr == chamber ) {
+                 //Momentum
+                 Double_t Px, Py, Pz, Pr;
+
+                 //Angle
+                 Double_t theta;
+                 Int_t    iTheta;
+                                                 
+                 Px = trackParam -> Px()   ;
+                 Py = trackParam -> Py()   ;
+                 Pz = trackParam -> Pz()   ;
+                 Pr = sqrt( Px*Px + Py*Py );
+
+                 //The calculation of theta, the angle of incidence:                                              
+                 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
+
+                 if      ( theta < 79 ) iTheta = 11;
+                 else if ( theta < 90 ) iTheta = int( theta - 79.);
+                 else                   iTheta = 11;
+
+                 eff    [iTheta] = eff    [iTheta] + 1;
+                 trackN [iTheta] = trackN [iTheta] + 1;
+                 chamber         = chamber + 1;
+               }
+
+               else {
+                 //Positions
+                 Double_t chamberZpos;
+
+                 //Momentum
+                 Double_t Px, Py, Pz, Pr;
+
+                 //Angles
+                 Double_t theta = 0.;
+                 Int_t    iTheta = 0 ;
+                                                  
+                 chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
+                 AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
+                                                  
+                 Px = trackParam -> Px()   ;
+                 Py = trackParam -> Py()   ;
+                 Pz = trackParam -> Pz()   ;
+                 Pr = sqrt( Px*Px + Py*Py );
+                                                 
+                 //The calculation of thetaI, the angle of incidence:
+                 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
+
+                 if      ( theta < 79 ) iTheta = 11;
+                 else if ( theta < 90 ) iTheta = int( theta - 79.);
+                 else                   iTheta = 11;
+                                                  
+                 eff    [iTheta] = eff    [iTheta] + 0;
+                 trackN [iTheta] = trackN [iTheta] + 1;
+                 chamber         = chamber + 1;
+                 iHit            = iHit - 1;
+               }
+             }
+           }
+         }
+                                                  
+         if ( iHit == nHits-1 ) {
+           if ( chamber == lastChamber )
+               continue;
+           else {
+             oldChamber = -1;
+             k          =  5;
+             iHit       = iHit-1;
+           }
+         }
+       }
+       //End loop on hits
+
+      } 
+      //End Loop on tracks
+
+      muondata.ResetRecTracks();
+      if ( event2Check!=0 )
+         iEvent=nEvents;
+    }
+    //End Loop on events
+
+    Double_t t [11];
+    Double_t efficiency [11];
+    Int_t i = 11;
+
+    Int_t th;
+    TGraph * effTheta = new TGraph ();
+
+    if ( firstChamber == lastChamber ) {
+      for ( th = 0; th < 11; ++th )
+      {
+       cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
+       cout<<"The chamber "<<firstChamber<<" has responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
+  
+       t [th]  = th + 79.5  ;
+       Double_t e = eff    [th];
+       Double_t n = trackN [th];
+       if ( n == 0. ) {
+         efficiency [th] = 0.;
+         cout<<"Efficiency = ? % (IS UNKNOWN)          \n";
+       }
+       else {
+         efficiency [th] = 100.*e/n;
+         cout<<"Efficiency = "<<efficiency [th]<<" %\n";
+       }
+      }
+    }
+
+    else{
+      for ( th = 0; th < 11; ++th )
+      {
+       cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
+       cout<<"The chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
+  
+       t [th]  = th + 79.5  ;
+       Double_t e = eff    [th];
+       Double_t n = trackN [th];
+       if ( n == 0. ) {
+         efficiency [th] = 0.;
+         cout<<"Efficiency = ? % (IS UNKNOWN)          \n";
+       }
+       else {
+         efficiency [th] = 100.*e/n;
+         cout<<"Efficiency = "<<efficiency [th]<<" %\n";
+       }
+      }
+    }
+
+    TCanvas * c1 = new TCanvas ();
+    effTheta = new TGraph( i, t, efficiency );
+
+    effTheta -> Draw("ALP");
+    MUONLoader->UnloadTracks();
+    c1->Update();
+}
+
+
+
+
+
+
+/*****************************************************************************************************************/
+/*****************************************************************************************************************/
+                                                /*RESOLUTION*/
+
+void resolution( Int_t event2Check=0, char * filename="galice.root" )
+{
+  cout<<"\nChamber(s) chosen;\nFirst chamber:";
+  cin>>firstChamber;
+  cout<<"Last chamber:";
+  cin>>lastChamber;
+  cout<<"\n\n";
+
+  TH1F * hDeltax;
+  TH1F * hDeltay;
+  TH2  * hDelta3D;
+
+  hDeltax    = new TH1F ("hDeltax " , "No Interest", 100, -0.3, 0.3);
+  hDeltay    = new TH1F ("hDeltay " , "No Interest", 500, -0.1, 0.1);
+  hDelta3D   = new TH2F ("hDelta3D" , "No Interest", 100, -0.3, 0.3, 500, -0.1, 0.1);
+
+
+//Creating Run Loader and openning file containing Hits
+//--------------------------------------------------------------------------------------------------------------
+  AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
+  if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
+
+
+//Getting MUONLoader
+//--------------------------------------------------------------------------------------------------------------
+  AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
+  MUONLoader -> LoadTracks("READ");
+  MUONLoader -> LoadRecPoints("READ");
+
+
+//Creating a MUON data container
+//--------------------------------------------------------------------------------------------------------------
+  AliMUONData muondata(MUONLoader,"MUON","MUON");
+
+
+  nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
+    
+    //Loop on events
+     for ( iEvent = 0; iEvent < nEvents; ++iEvent )
+    { 
+      if (event2Check!=0)
+         iEvent=event2Check;
+      printf("\r>>> Event %d ",iEvent);
+      RunLoader->GetEvent(iEvent);
+      
+      //Addressing
+      muondata.SetTreeAddress("RT");     
+      muondata.GetRecTracks();
+      tracks = muondata.RecTracks();
+
+      //Loop on track
+      nTracks = (Int_t) tracks -> GetEntriesFast();
+      for ( iTrack = 0; iTrack < nTracks; ++iTrack )
+      {
+       track       = (AliMUONTrack*) tracks ->At(iTrack)          ;
+       trackParams =                 track  ->GetTrackParamAtHit();
+       hits        =                 track  ->GetHitForRecAtHit() ;
+       
+       //Loop on hits
+       nHits = (Int_t) hits -> GetEntriesFast();
+       for ( iHit = 0; iHit < nHits; ++iHit )
+       { 
+         trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
+         hit        = (AliMUONHitForRec*)  hits        -> At(iHit);
+         chamberNbr = hit -> GetChamberNumber();
+         if ( chamberNbr >= firstChamber-1 ) {
+           if ( chamberNbr < lastChamber ) {
+             //Positions
+             Double_t traX, traY;
+             Double_t hitX, hitY;
+                                                    
+             //Resolution
+             Double_t deltaX, deltaY;
+                                                    
+             traX = trackParam -> GetNonBendingCoor();
+             traY = trackParam -> GetBendingCoor();
+             hitX = hit        -> GetNonBendingCoor();
+             hitY = hit        -> GetBendingCoor();
+                                                    
+             deltaX = traX - hitX;
+             deltaY = traY - hitY;
+                                                    
+             hDeltax  -> Fill (deltaX);
+             hDeltay  -> Fill (deltaY);
+             hDelta3D -> Fill (deltaX, deltaY);
+           }
+         }
+       }
+       //End loop on hits
+      }
+      //End loop on tracks
+      muondata.ResetRecTracks();
+      if ( event2Check!=0 )
+         iEvent=nEvents;
+    }
+    //End loop on events
+
+    char hXtitle[80];
+    char hYtitle[80];
+    char h3title[80];
+
+    if ( firstChamber == lastChamber ) {
+      sprintf(hXtitle,"DeltaX Chamber %d",firstChamber);
+      sprintf(hYtitle,"DeltaY Chamber %d",lastChamber);
+      sprintf(h3title,"DeltaX and Delta Y Chamber %d",lastChamber);
+    }
+    else{
+      sprintf(hXtitle,"DeltaX Chamber %d to %d",firstChamber,lastChamber);
+      sprintf(hYtitle,"DeltaY Chamber %d to %d",firstChamber,lastChamber);
+      sprintf(h3title,"DeltaX and Delta Y Chamber %d to %d",firstChamber,lastChamber);
+    }
+  
+
+    hDeltax  -> SetTitle (hXtitle);
+    hDeltay  -> SetTitle (hYtitle);
+    hDelta3D -> SetTitle (h3title);
+
+    Double_t rmsX = hDeltax -> GetRMS     ();
+    Double_t errX = hDeltax -> GetRMSError();
+
+    TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
+    hDeltay -> Fit("fY","R","E");
+    Double_t sigY = fY -> GetParameter(2);
+    Double_t errY = fY -> GetParError (2);
+
+    if ( firstChamber == lastChamber ) {
+      cout<<"\nThe resolution (Non Bending plane) of the chamber "<<firstChamber<<" = "<<rmsX<<" cm +-"<<errX;
+      cout<<"\nThe resolution (Bending plane) of the chamber "<<firstChamber<<" = "<<sigY<<" cm +- "<<errY;
+    }
+    else {
+      cout<<"\nThe resolution (Non Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<rmsX<<" cm +- "<<errX;
+      cout<<"\nThe resolution (Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<sigY<<" cm +- "<<errY;
+    }
+
+    cout<<"\n\n\n";
+    MUONLoader->UnloadTracks();
+
+}
+
+
+
+
+
+
+
+/*****************************************************************************************************************/
+/*****************************************************************************************************************/
+                                       /*RESOLUTION IN FUNCTION OF PHI*/
+
+void resolutionPhi( Int_t event2Check=0, char * filename="galice.root" )
+{
+  cout<<"\nChamber(s) chosen;\nFirst chamber:";
+  cin>>firstChamber;
+  cout<<"Last chamber:";
+  cin>>lastChamber;
+  cout<<"\n\n";
+
+  TH1F * hDeltaX[5];
+  TH1F * hDeltaY[5];
+
+  hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
+  
+  hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
+
+
+//Creating Run Loader and openning file containing Hits
+//--------------------------------------------------------------------------------------------------------------
+  AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
+  if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
+
+
+//Getting MUONLoader
+//--------------------------------------------------------------------------------------------------------------
+  AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
+  MUONLoader -> LoadTracks("READ");
+  MUONLoader -> LoadRecPoints("READ");
+
+
+//Creating a MUON data container
+//--------------------------------------------------------------------------------------------------------------
+  AliMUONData muondata(MUONLoader,"MUON","MUON");
+
+
+  nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
+
+    //Loop on events
+    for ( iEvent = 0; iEvent < nEvents; ++iEvent )
+    {
+      if ( event2Check!=0 )
+         iEvent=event2Check;
+      printf("\r>>> Event %d ",iEvent);
+      RunLoader->GetEvent(iEvent);
+      
+      //Addressing
+      muondata.SetTreeAddress("RT");     
+      muondata.GetRecTracks();
+      tracks = muondata.RecTracks();
+
+      //Loop on track
+      nTracks = (Int_t) tracks -> GetEntriesFast();
+      for ( iTrack = 0; iTrack < nTracks; ++iTrack )
+      {
+       track       = (AliMUONTrack*) tracks ->At(iTrack)          ;
+       trackParams =                 track  ->GetTrackParamAtHit();
+       hits        =                 track  ->GetHitForRecAtHit() ;
+       
+       //Loop on hits
+       nHits = (Int_t) hits -> GetEntriesFast();
+       for ( iHit = 0; iHit < nHits; ++iHit )
+       { 
+         trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
+         hit        = (AliMUONHitForRec* ) hits        -> At(iHit);
+         chamberNbr = hit -> GetChamberNumber();
+         if ( chamberNbr >= firstChamber -1 ) {
+           if ( chamberNbr < lastChamber     ) {
+             //Positions
+             Double_t traX, traY;
+             Double_t hitX, hitY;
+             Double_t aveY, aveX;
+
+             //Angles
+             Double_t phi;
+             Int_t    iPhi;
+
+             //Resolution                                                                           
+             Double_t deltaX;
+             Double_t deltaY;
+
+             traX = trackParam -> GetNonBendingCoor();
+             traY = trackParam -> GetBendingCoor()   ;
+                                                                                                     
+             hitX = hit        -> GetNonBendingCoor();
+             hitY = hit        -> GetBendingCoor()   ;
+                                                                                                     
+             aveX = 1./2.*(traX + hitX);
+             aveY = 1./2.*(traY + hitY);
+
+             //The calculation of phi: 
+             phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));                                             
+
+             if ( phi < 0.) phi = 360 - abs(phi);
+             iPhi = int( phi/72. );
+                                                     
+             deltaX = traX - hitX;
+             deltaY = traY - hitY;
+
+             hDeltaX [iPhi] -> Fill( deltaX );
+             hDeltaY [iPhi] -> Fill( deltaY );
+           }
+         }
+           
+       }
+       //End loop on hits
+
+      }
+      //End loop on tracks
+      muondata.ResetRecTracks();
+      if ( event2Check!=0 )
+         iEvent=nEvents;
+
+    }
+    //End loop on events
+
+
+    Int_t iPhi;
+    Int_t iPhiMax = 5;
+    Int_t phiMin, phiMax;
+    Float_t phi[5];
+    Float_t sigmaY[5];
+    Float_t errSY [5];
+    Float_t rmsX  [5];
+    Float_t errSX [5];
+    Float_t errPh [5];
+
+    for ( iPhi = 0; iPhi < iPhiMax; ++iPhi )
+    {
+      char hXtitle[80];
+      char hYtitle[80];
+
+      phiMin = 72*iPhi     ;
+      phiMax = 72*iPhi + 72;
+        
+      TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
+       
+      if ( firstChamber == lastChamber ) {
+       sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d",phiMin,phiMax,firstChamber);
+       sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d",phiMin,phiMax,lastChamber);
+      }
+      else{
+       sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
+       sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
+      }
+
+      hDeltaY [iPhi] -> SetTitle(hYtitle);
+      hDeltaX [iPhi] -> SetTitle(hXtitle);
+
+      hDeltaY [iPhi]      -> Fit("fY","R","E");
+      sigmaY  [iPhi] = fY -> GetParameter(2)  ;
+      errSY   [iPhi] = fY -> GetParError(2)   ;
+
+      rmsX    [iPhi] = hDeltaX [iPhi] -> GetRMS();
+      errSX   [iPhi] = hDeltaX [iPhi] -> GetRMSError(); 
+
+      phi     [iPhi] = 72*iPhi + 36 ;
+      errPh   [iPhi] = 36;
+    }
+
+    //--------------------------------------------------------------------------------------------------------------
+    //For plotting resolution in function of the angle of incidence
+
+    TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
+    c1-> Divide(1,2);
+    c1->cd(1);
+
+    TGraphErrors * GraphX = new TGraphErrors( iPhiMax, phi, rmsX, errPh, errSX);
+    GraphX->SetTitle("Resolution en X (Phi)");
+    GraphX->Draw("ALP");
+
+    c1->cd(2);
+    TGraphErrors * GraphY = new TGraphErrors( iPhiMax, phi, sigmaY, errPh, errSY);
+    GraphY->SetTitle("Resolution en Y (Phi)");
+    GraphY->Draw("ALP");
+
+    cout<<"\n\n\n";
+
+    MUONLoader->UnloadTracks();
+
+}
+
+
+
+
+
+
+
+/*****************************************************************************************************************/
+/*****************************************************************************************************************/
+                                       /*RESOLUTION IN FUNCTION OF THETA*/
+
+void resolutionTheta( Int_t event2Check=0, char * filename="galice.root" )
+{
+  cout<<"\nChamber(s) chosen;\nFirst chamber:";
+  cin>>firstChamber;
+  cout<<"Last chamber:";
+  cin>>lastChamber;
+  cout<<"\n\n";
+
+  TH1F * hDeltaX[9];
+  TH1F * hDeltaY[9];
+
+  hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
+  
+  hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
+
+
+//Creating Run Loader and openning file containing Hits
+//--------------------------------------------------------------------------------------------------------------
+  AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
+  if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
+
+
+//Getting MUONLoader
+//--------------------------------------------------------------------------------------------------------------
+  AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
+  MUONLoader -> LoadTracks("READ");
+  MUONLoader -> LoadRecPoints("READ");
+
+
+//Creating a MUON data container
+//--------------------------------------------------------------------------------------------------------------
+  AliMUONData muondata(MUONLoader,"MUON","MUON");
+
+
+  nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
+
+    //Loop on events
+    for ( iEvent = 0; iEvent < nEvents; ++iEvent )
+    {
+      if ( event2Check!=0 )
+         iEvent=event2Check;
+      printf("\r>>> Event %d ",iEvent);
+      RunLoader->GetEvent(iEvent);
+      
+      //Addressing
+      muondata.SetTreeAddress("RT");     
+      muondata.GetRecTracks();
+      tracks = muondata.RecTracks();
+
+      //Loop on track
+      nTracks = (Int_t) tracks -> GetEntriesFast();
+      for ( iTrack = 0; iTrack < nTracks; ++iTrack )
+      { 
+       track       = (AliMUONTrack*) tracks ->At(iTrack)          ;
+       trackParams =                 track  ->GetTrackParamAtHit();
+       hits        =                 track  ->GetHitForRecAtHit() ;
+       
+       //Loop on hits
+       nHits = (Int_t) hits -> GetEntriesFast();
+       for ( iHit = 0; iHit < nHits; ++iHit )
+       { 
+         trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
+         hit        = (AliMUONHitForRec* ) hits        -> At(iHit);
+         chamberNbr = hit -> GetChamberNumber();
+         if ( chamberNbr >= firstChamber -1 ) {
+           if ( chamberNbr < lastChamber ) {
+             //Positions
+             Double_t traX, traY;
+             Double_t hitX, hitY, hitZ;
+             Double_t aveY, aveX;
+                                                     
+             //Angles
+             Double_t theta;
+             Int_t    iTheta;
+
+             //Resolution                                                                           
+             Double_t deltaX;
+             Double_t deltaY;
+
+             traX = trackParam -> GetNonBendingCoor();
+             traY = trackParam -> GetBendingCoor()   ;
+                                                                                                     
+             hitX = hit        -> GetNonBendingCoor();
+             hitY = hit        -> GetBendingCoor()   ;
+             hitZ = hit        -> GetZ();
+                                                     
+             aveX = 1./2.*(traX + hitX);
+             aveY = 1./2.*(traY + hitY);
+                                                     
+             //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
+             theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ ));
+
+             iTheta = int( theta );
+             if( theta > 10 ) iTheta = 9;
+             if( theta < 1  ) iTheta = 1; 
+                                                     
+             deltaX = traX - hitX;
+             deltaY = traY - hitY;
+
+             hDeltaX [iTheta-1] -> Fill( deltaX );
+             hDeltaY [iTheta-1] -> Fill( deltaY );
+           }
+         }
+                                                   
+       }
+       //End loop on hits
+
+      }
+      //End loop on tracks
+      muondata.ResetRecTracks();
+      if ( event2Check!=0 )
+         iEvent=nEvents;
+
+    }
+    //End loop on events
+
+
+    Int_t iTheta;
+    Int_t iThetaMax = 9;
+    Int_t thetaMin, thetaMax;
+    Float_t theta [9];
+    Float_t sigmaY[9];
+    Float_t errSY [9];
+    Float_t rmsX  [9];
+    Float_t errSX [9];
+    Float_t ErrTh [9];
+
+    for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
+    {
+      char hXtitle[80];
+      char hYtitle[80];
+
+      //To find the polar angle
+      thetaMin = 178 - iTheta;
+      thetaMax = 179 - iTheta;
+        
+      TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
+       
+      if ( firstChamber == lastChamber ) {
+       sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
+       sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
+      }
+      else{
+       sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
+       sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
+      }
+
+      hDeltaY [iTheta] -> SetTitle(hYtitle);
+      hDeltaX [iTheta] -> SetTitle(hXtitle);
+
+      hDeltaY [iTheta]      -> Fit("fY","R","E");
+      sigmaY  [iTheta] = fY -> GetParameter(2)  ;
+      errSY   [iTheta] = fY -> GetParError(2)   ;
+
+      rmsX    [iTheta] = hDeltaX [iTheta] -> GetRMS();
+      errSX   [iTheta] = hDeltaX [iTheta] -> GetRMSError(); 
+
+      theta   [iTheta] = 178.5 - iTheta;
+      ErrTh   [iTheta] = 0.5;
+    }
+
+    //--------------------------------------------------------------------------------------------------------------
+    //For plotting resolution in function of the angle of incidence
+
+    TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
+    c1-> Divide(1,2);
+    c1->cd(1);
+
+    TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, ErrTh, errSX);
+    GraphX->SetTitle("Resolution en X (Theta)");
+    GraphX->Draw("ALP");
+
+    c1->cd(2);
+    TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, ErrTh, errSY);
+    GraphY->SetTitle("Resolution en Y (Theta)");
+    GraphY->Draw("ALP");
+
+    cout<<"\n\n\n";
+    MUONLoader->UnloadTracks();
+
+}
+
+
+
+
+
+/*****************************************************************************************************************/
+/*****************************************************************************************************************/
+                              /*RESOLUTION IN FUNCTION OF THETA OF INCIDENCE*/
+
+void resolutionThetaI( Int_t event2Check=0, char * filename="galice.root" )
+{
+  cout<<"\nChamber(s) chosen;\nFirst chamber:";
+  cin>>firstChamber;
+  cout<<"Last chamber:";
+  cin>>lastChamber;
+  cout<<"\n\n";
+
+  TH1F * hDeltaX[11];
+  TH1F * hDeltaY[11];
+
+  hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[9] = new TH1F ("hDeltaX9" , "No Interest", 100, -0.3, 0.3);
+  hDeltaX[10]= new TH1F ("hDeltaX10", "No Interest", 100, -0.3, 0.3);
+  
+  hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[9] = new TH1F ("hDeltaY9" , "No Interest", 500, -0.1, 0.1);
+  hDeltaY[10]= new TH1F ("hDeltaY10", "No Interest", 500, -0.1, 0.1);
+
+
+//Creating Run Loader and openning file containing Hits
+//--------------------------------------------------------------------------------------------------------------
+  AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
+  if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
+
+
+//Getting MUONLoader
+//--------------------------------------------------------------------------------------------------------------
+  AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
+  MUONLoader -> LoadTracks("READ");
+  MUONLoader -> LoadRecPoints("READ");
+
+
+//Creating a MUON data container
+//--------------------------------------------------------------------------------------------------------------
+  AliMUONData muondata(MUONLoader,"MUON","MUON");
+
+
+  nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
+
+    //Loop on events
+    for ( iEvent = 0; iEvent < nEvents; ++iEvent )
+    { 
+      if ( event2Check!=0 )
+         iEvent=event2Check;
+      printf("\r>>> Event %d ",iEvent);
+      RunLoader->GetEvent(iEvent);
+      
+      //Addressing
+      muondata.SetTreeAddress("RT");     
+      muondata.GetRecTracks();
+      tracks = muondata.RecTracks();
+
+      //Loop on track
+      nTracks = (Int_t) tracks -> GetEntriesFast();
+      for ( iTrack = 0; iTrack < nTracks; ++iTrack )
+      {
+       track       = (AliMUONTrack*) tracks ->At(iTrack)          ;
+       trackParams =                 track  ->GetTrackParamAtHit();
+       hits        =                 track  ->GetHitForRecAtHit() ;
+       
+       //Loop on hits
+       nHits = (Int_t) hits -> GetEntriesFast();
+       for ( iHit = 0; iHit < nHits; ++iHit )
+       {
+         trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
+         hit        = (AliMUONHitForRec* ) hits        -> At(iHit);
+         chamberNbr = hit -> GetChamberNumber();
+         if ( chamberNbr >= firstChamber - 1 ) {
+           if ( chamberNbr < lastChamber ) {
+             //Momentum
+             Double_t Px, Py, Pz, Pr;
+
+             //Positions
+             Double_t traX, traY;
+             Double_t hitX, hitY;
+
+             //Resolution
+             Double_t deltaX;
+             Double_t deltaY;
+
+             //Angle
+             Double_t theta;
+             Int_t    iTheta;
+
+             Px = trackParam -> Px()   ;
+             Py = trackParam -> Py()   ;
+             Pz = trackParam -> Pz()   ;
+             Pr = sqrt( Px*Px + Py*Py );
+
+             traX = trackParam -> GetNonBendingCoor();
+             traY = trackParam -> GetBendingCoor();
+
+             hitX = hit        -> GetNonBendingCoor();
+             hitY = hit        -> GetBendingCoor();                                        
+
+             //The calculation of theta, the angle of incidence:
+             theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr))
+;
+             if ( theta < 79 ) iTheta = 0;
+             else iTheta = int( theta - 79.);
+
+             deltaX = traX - hitX;
+             deltaY = traY - hitY;
+
+             hDeltaX [iTheta] -> Fill (deltaX);
+             hDeltaY [iTheta] -> Fill (deltaY);
+           }
+         }
+                                                   
+       }
+       //End loop on hits
+
+      }
+      //End loop on tracks
+      muondata.ResetRecTracks();
+      if ( event2Check!=0 )
+         iEvent=nEvents;
+
+    }
+    //End loop on events
+
+
+    Int_t iTheta;
+    Int_t iThetaMax = 11;
+    Int_t thetaMin, thetaMax;
+    Float_t theta [11];
+    Float_t sigmaY[11];
+    Float_t errSY [11];
+    Float_t rmsX  [11];
+    Float_t errSX [11];
+    Float_t errTh [11];
+
+    for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
+    {
+      char hXtitle[80];
+      char hYtitle[80];
+
+      thetaMin = iTheta + 79;
+      thetaMax = iTheta + 80;
+        
+      TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
+       
+      if ( firstChamber == lastChamber ) {
+       sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
+       sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
+      }
+      else{
+       sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
+       sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
+      }
+
+      hDeltaY [iTheta] -> SetTitle(hYtitle);
+      hDeltaX [iTheta] -> SetTitle(hXtitle);
+
+      hDeltaY [iTheta]      -> Fit("fY","R","E");
+      sigmaY  [iTheta] = fY -> GetParameter(2)  ;
+      errSY   [iTheta] = fY -> GetParError(2)   ;
+
+      rmsX    [iTheta] = hDeltaX [iTheta] -> GetRMS();
+      errSX   [iTheta] = hDeltaX [iTheta] -> GetRMSError(); 
+
+      theta   [iTheta] = iTheta + 79.5 ;
+      errTh   [iTheta] = 0.5;
+    }
+
+    //--------------------------------------------------------------------------------------------------------------
+    //For plotting resolution in function of the angle of incidence
+
+    TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
+    c1-> Divide(1,2);
+    c1->cd(1);
+
+    TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, errTh, errSX);
+    GraphX->SetTitle("Resolution en X (Theta)");
+    GraphX->Draw("ALP");
+
+    c1->cd(2);
+    TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, errTh, errSY);
+    GraphY->SetTitle("Resolution en Y (Theta)");
+    GraphY->Draw("ALP");
+
+    cout<<"\n\n\n";
+    MUONLoader->UnloadTracks();
+
+}
+
+
+
+