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 "AliMagFMaps.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   printf("Loading field map...\n");
 249   if (!AliTracker::GetFieldMap()) {
 250   AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
 251   AliTracker::SetFieldMap(field, kFALSE);}
 252
 253
 254 //--------------------------------------------------------------------------------------------------------------
 255 //Set Field Map for track extrapolation
 256   AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
 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   printf("Loading field map...\n");
 514   if (!AliTracker::GetFieldMap()) {
 515   AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
 516   AliTracker::SetFieldMap(field, kFALSE);}
 517
 518
 519 //--------------------------------------------------------------------------------------------------------------
 520 //Set Field Map for track extrapolation
 521   AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
 522
 523
 524 //Creating a MUON data container
 525 //--------------------------------------------------------------------------------------------------------------
 526   AliMUONData muondata(MUONLoader,"MUON","MUON");
 527
 528
 529   nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
 530
 531     //Loop on events
 532     for ( iEvent = 0; iEvent < nEvents; ++iEvent )
 533     {
 534       if ( event2Check!=0 )
 535           iEvent=event2Check;
 536       printf("\r>>> Event %d ",iEvent);
 537       RunLoader->GetEvent(iEvent);
 538
 539       //Addressing
 540       muondata.SetTreeAddress("RT");
 541       muondata.GetRecTracks();
 542       tracks = muondata.RecTracks();
 543
 544       //Loop on track
 545       nTracks = (Int_t) tracks -> GetEntriesFast();
 546       for ( iTrack = 0; iTrack < nTracks; ++iTrack )
 547       {
 548         track       = (AliMUONTrack*) tracks ->At(iTrack)          ;
 549         trackParams =                 track  ->GetTrackParamAtCluster();
 550         hits        =                 track  ->GetHitForRecAtHit() ;
 551         chamber     = firstChamber - 1;
 552         oldChamber  = -1;
 553         Int_t k     =  1;
 554
 555         //Loop on hits
 556         nHits = (Int_t) hits -> GetEntriesFast();
 557         for ( iHit = 0; iHit < nHits; ++iHit )
 558         {
 559           trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
 560           hit        = (AliMUONHitForRec*)  hits        -> At(iHit);
 561           chamberNbr = hit -> GetChamberNumber();
 562
 563           if ( chamberNbr == oldChamber )
 564               continue;
 565           else {
 566             oldChamber = chamberNbr;
 567             if ( chamberNbr > chamber - k ) {
 568               if ( chamber < lastChamber ) {
 569                 if ( chamberNbr == chamber ) {
 570                   //Momentum
 571                   Double_t Px, Py, Pz, Pr;
 572
 573                   //Angle
 574                   Double_t theta;
 575                   Int_t    iTheta;
 576
 577                   Px = trackParam -> Px()   ;
 578                   Py = trackParam -> Py()   ;
 579                   Pz = trackParam -> Pz()   ;
 580                   Pr = sqrt( Px*Px + Py*Py );
 581
 582                   //The calculation of theta, the angle of incidence:
 583                   theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
 584
 585                   if      ( theta < 79 ) iTheta = 11;
 586                   else if ( theta < 90 ) iTheta = int( theta - 79.);
 587                   else                   iTheta = 11;
 588
 589                   eff    [iTheta] = eff    [iTheta] + 1;
 590                   trackN [iTheta] = trackN [iTheta] + 1;
 591                   chamber         = chamber + 1;
 592                 }
 593
 594                 else {
 595                   //Positions
 596                   Double_t chamberZpos;
 597
 598                   //Momentum
 599                   Double_t Px, Py, Pz, Pr;
 600
 601                   //Angles
 602                   Double_t theta = 0.;
 603                   Int_t    iTheta = 0 ;
 604
 605                   chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
 606                   AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
 607
 608                   Px = trackParam -> Px()   ;
 609                   Py = trackParam -> Py()   ;
 610                   Pz = trackParam -> Pz()   ;
 611                   Pr = sqrt( Px*Px + Py*Py );
 612
 613                   //The calculation of thetaI, the angle of incidence:
 614                   theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
 615
 616                   if      ( theta < 79 ) iTheta = 11;
 617                   else if ( theta < 90 ) iTheta = int( theta - 79.);
 618                   else                   iTheta = 11;
 619
 620                   eff    [iTheta] = eff    [iTheta] + 0;
 621                   trackN [iTheta] = trackN [iTheta] + 1;
 622                   chamber         = chamber + 1;
 623                   iHit            = iHit - 1;
 624                 }
 625               }
 626             }
 627           }
 628
 629           if ( iHit == nHits-1 ) {
 630             if ( chamber == lastChamber )
 631                 continue;
 632             else {
 633               oldChamber = -1;
 634               k          =  5;
 635               iHit       = iHit-1;
 636             }
 637           }
 638         }
 639         //End loop on hits
 640
 641       }
 642       //End Loop on tracks
 643
 644       muondata.ResetRecTracks();
 645       if ( event2Check!=0 )
 646           iEvent=nEvents;
 647     }
 648     //End Loop on events
 649
 650     Double_t t [11];
 651     Double_t efficiency [11];
 652     Int_t i = 11;
 653
 654     Int_t th;
 655     TGraph * effTheta = new TGraph ();
 656
 657     if ( firstChamber == lastChamber ) {
 658       for ( th = 0; th < 11; ++th )
 659       {
 660         cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
 661         cout<<"The chamber "<<firstChamber<<" has responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
 662
 663         t [th]  = th + 79.5  ;
 664         Double_t e = eff    [th];
 665         Double_t n = trackN [th];
 666
 667         if ( n == 0. ) {
 668           efficiency [th] = 0.;
 669           cout<<"Efficiency = ? % (IS UNKNOWN)          \n";
 670         }
 671         else {
 672           efficiency [th] = 100.*e/n;
 673           cout<<"Efficiency = "<<efficiency [th]<<" %\n";
 674         }
 675       }
 676     }
 677
 678     else{
 679       for ( th = 0; th < 11; ++th )
 680       {
 681         cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
 682         cout<<"The chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
 683
 684         t [th]  = th + 79.5  ;
 685         Double_t e = eff    [th];
 686         Double_t n = trackN [th];
 687
 688         if ( n == 0. ) {
 689           efficiency [th] = 0.;
 690           cout<<"Efficiency = ? % (IS UNKNOWN)          \n";
 691         }
 692         else {
 693           efficiency [th] = 100.*e/n;
 694           cout<<"Efficiency = "<<efficiency [th]<<" %\n";
 695         }
 696       }
 697     }
 698
 699     TCanvas * c1 = new TCanvas ();
 700     effTheta = new TGraph( i, t, efficiency );
 701
 702     effTheta -> Draw("ALP");
 703
 704     MUONLoader->UnloadTracks();
 705     c1->Update();
 706 }
 707
 708
 709
 710
 711
 712
 713 /*****************************************************************************************************************/
 714 /*****************************************************************************************************************/
 715                                                 /*RESOLUTION*/
 716
 717 void resolution( Int_t event2Check=0, char * filename="galice.root" )
 718 {
 719   cout<<"\nChamber(s) chosen;\nFirst chamber:";
 720   cin>>firstChamber;
 721   cout<<"Last chamber:";
 722   cin>>lastChamber;
 723   cout<<"\n\n";
 724
 725   TH1F * hDeltax;
 726   TH1F * hDeltay;
 727   TH2  * hDelta3D;
 728
 729   hDeltax    = new TH1F ("hDeltax " , "No Interest", 100, -0.3, 0.3);
 730   hDeltay    = new TH1F ("hDeltay " , "No Interest", 500, -0.1, 0.1);
 731   hDelta3D   = new TH2F ("hDelta3D" , "No Interest", 100, -0.3, 0.3, 500, -0.1, 0.1);
 732
 733
 734 //Creating Run Loader and openning file containing Hits
 735 //--------------------------------------------------------------------------------------------------------------
 736   AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
 737   if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
 738
 739
 740 //Getting MUONLoader
 741 //--------------------------------------------------------------------------------------------------------------
 742   AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
 743   MUONLoader -> LoadTracks("READ");
 744   MUONLoader -> LoadRecPoints("READ");
 745
 746
 747 //Creating a MUON data container
 748 //--------------------------------------------------------------------------------------------------------------
 749   AliMUONData muondata(MUONLoader,"MUON","MUON");
 750
 751
 752   nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
 753
 754     //Loop on events
 755      for ( iEvent = 0; iEvent < nEvents; ++iEvent )
 756     {
 757       if (event2Check!=0)
 758           iEvent=event2Check;
 759       printf("\r>>> Event %d ",iEvent);
 760       RunLoader->GetEvent(iEvent);
 761
 762       //Addressing
 763       muondata.SetTreeAddress("RT");
 764       muondata.GetRecTracks();
 765       tracks = muondata.RecTracks();
 766
 767       //Loop on track
 768       nTracks = (Int_t) tracks -> GetEntriesFast();
 769       for ( iTrack = 0; iTrack < nTracks; ++iTrack )
 770       {
 771         track       = (AliMUONTrack*) tracks ->At(iTrack)          ;
 772         trackParams =                 track  ->GetTrackParamAtCluster();
 773         hits        =                 track  ->GetHitForRecAtHit() ;
 774
 775         //Loop on hits
 776         nHits = (Int_t) hits -> GetEntriesFast();
 777         for ( iHit = 0; iHit < nHits; ++iHit )
 778         {
 779           trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
 780           hit        = (AliMUONHitForRec*)  hits        -> At(iHit);
 781           chamberNbr = hit -> GetChamberNumber();
 782           if ( chamberNbr >= firstChamber-1 ) {
 783             if ( chamberNbr < lastChamber ) {
 784               //Positions
 785               Double_t traX, traY;
 786               Double_t hitX, hitY;
 787
 788               //Resolution
 789               Double_t deltaX, deltaY;
 790
 791               traX = trackParam -> GetNonBendingCoor();
 792               traY = trackParam -> GetBendingCoor();
 793               hitX = hit        -> GetNonBendingCoor();
 794               hitY = hit        -> GetBendingCoor();
 795
 796               deltaX = traX - hitX;
 797               deltaY = traY - hitY;
 798
 799               hDeltax  -> Fill (deltaX);
 800               hDeltay  -> Fill (deltaY);
 801               hDelta3D -> Fill (deltaX, deltaY);
 802             }
 803           }
 804         }
 805         //End loop on hits
 806       }
 807       //End loop on tracks
 808       muondata.ResetRecTracks();
 809       if ( event2Check!=0 )
 810           iEvent=nEvents;
 811     }
 812     //End loop on events
 813
 814     char hXtitle[80];
 815     char hYtitle[80];
 816     char h3title[80];
 817
 818     if ( firstChamber == lastChamber ) {
 819       sprintf(hXtitle,"DeltaX Chamber %d",firstChamber);
 820       sprintf(hYtitle,"DeltaY Chamber %d",lastChamber);
 821       sprintf(h3title,"DeltaX and Delta Y Chamber %d",lastChamber);
 822     }
 823     else{
 824       sprintf(hXtitle,"DeltaX Chamber %d to %d",firstChamber,lastChamber);
 825       sprintf(hYtitle,"DeltaY Chamber %d to %d",firstChamber,lastChamber);
 826       sprintf(h3title,"DeltaX and Delta Y Chamber %d to %d",firstChamber,lastChamber);
 827     }
 828
 829
 830     hDeltax  -> SetTitle (hXtitle);
 831     hDeltay  -> SetTitle (hYtitle);
 832     hDelta3D -> SetTitle (h3title);
 833
 834     Double_t rmsX = hDeltax -> GetRMS     ();
 835     Double_t errX = hDeltax -> GetRMSError();
 836
 837     TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
 838     hDeltay -> Fit("fY","R","E");
 839     Double_t sigY = fY -> GetParameter(2);
 840     Double_t errY = fY -> GetParError (2);
 841
 842     if ( firstChamber == lastChamber ) {
 843       cout<<"\nThe resolution (Non Bending plane) of the chamber "<<firstChamber<<" = "<<rmsX<<" cm +-"<<errX;
 844       cout<<"\nThe resolution (Bending plane) of the chamber "<<firstChamber<<" = "<<sigY<<" cm +- "<<errY;
 845     }
 846     else {
 847       cout<<"\nThe resolution (Non Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<rmsX<<" cm +- "<<errX;
 848       cout<<"\nThe resolution (Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<sigY<<" cm +- "<<errY;
 849     }
 850
 851     cout<<"\n\n\n";
 852     MUONLoader->UnloadTracks();
 853
 854 }
 855
 856
 857
 858
 859
 860
 861
 862 /*****************************************************************************************************************/
 863 /*****************************************************************************************************************/
 864                                        /*RESOLUTION IN FUNCTION OF PHI*/
 865
 866 void resolutionPhi( Int_t event2Check=0, char * filename="galice.root" )
 867 {
 868   cout<<"\nChamber(s) chosen;\nFirst chamber:";
 869   cin>>firstChamber;
 870   cout<<"Last chamber:";
 871   cin>>lastChamber;
 872   cout<<"\n\n";
 873
 874   TH1F * hDeltaX[5];
 875   TH1F * hDeltaY[5];
 876
 877   hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
 878   hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
 879   hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
 880   hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
 881   hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
 882
 883   hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
 884   hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
 885   hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
 886   hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
 887   hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
 888
 889
 890 //Creating Run Loader and openning file containing Hits
 891 //--------------------------------------------------------------------------------------------------------------
 892   AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
 893   if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
 894
 895
 896 //Getting MUONLoader
 897 //--------------------------------------------------------------------------------------------------------------
 898   AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
 899   MUONLoader -> LoadTracks("READ");
 900   MUONLoader -> LoadRecPoints("READ");
 901
 902
 903 //Creating a MUON data container
 904 //--------------------------------------------------------------------------------------------------------------
 905   AliMUONData muondata(MUONLoader,"MUON","MUON");
 906
 907
 908   nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
 909
 910     //Loop on events
 911     for ( iEvent = 0; iEvent < nEvents; ++iEvent )
 912     {
 913       if ( event2Check!=0 )
 914           iEvent=event2Check;
 915       printf("\r>>> Event %d ",iEvent);
 916       RunLoader->GetEvent(iEvent);
 917
 918       //Addressing
 919       muondata.SetTreeAddress("RT");
 920       muondata.GetRecTracks();
 921       tracks = muondata.RecTracks();
 922
 923       //Loop on track
 924       nTracks = (Int_t) tracks -> GetEntriesFast();
 925       for ( iTrack = 0; iTrack < nTracks; ++iTrack )
 926       {
 927         track       = (AliMUONTrack*) tracks ->At(iTrack)          ;
 928         trackParams =                 track  ->GetTrackParamAtCluster();
 929         hits        =                 track  ->GetHitForRecAtHit() ;
 930
 931         //Loop on hits
 932         nHits = (Int_t) hits -> GetEntriesFast();
 933         for ( iHit = 0; iHit < nHits; ++iHit )
 934         {
 935           trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
 936           hit        = (AliMUONHitForRec* ) hits        -> At(iHit);
 937           chamberNbr = hit -> GetChamberNumber();
 938           if ( chamberNbr >= firstChamber -1 ) {
 939             if ( chamberNbr < lastChamber     ) {
 940               //Positions
 941               Double_t traX, traY;
 942               Double_t hitX, hitY;
 943               Double_t aveY, aveX;
 944
 945               //Angles
 946               Double_t phi;
 947               Int_t    iPhi;
 948
 949               //Resolution
 950               Double_t deltaX;
 951               Double_t deltaY;
 952
 953               traX = trackParam -> GetNonBendingCoor();
 954               traY = trackParam -> GetBendingCoor()   ;
 955
 956               hitX = hit        -> GetNonBendingCoor();
 957               hitY = hit        -> GetBendingCoor()   ;
 958
 959               aveX = 1./2.*(traX + hitX);
 960               aveY = 1./2.*(traY + hitY);
 961
 962               //The calculation of phi:
 963               phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));
 964
 965               if ( phi < 0.) phi = 360. - fabs(phi);
 966               iPhi = int( phi/72. );
 967
 968               deltaX = traX - hitX;
 969               deltaY = traY - hitY;
 970
 971               hDeltaX [iPhi] -> Fill( deltaX );
 972               hDeltaY [iPhi] -> Fill( deltaY );
 973             }
 974           }
 975
 976         }
 977         //End loop on hits
 978
 979       }
 980       //End loop on tracks
 981       muondata.ResetRecTracks();
 982       if ( event2Check!=0 )
 983           iEvent=nEvents;
 984
 985     }
 986     //End loop on events
 987
 988
 989     Int_t iPhi;
 990     Int_t iPhiMax = 5;
 991     Int_t phiMin, phiMax;
 992
 993     Float_t phi[5];
 994     Float_t sigmaY[5];
 995     Float_t errSY [5];
 996     Float_t rmsX  [5];
 997     Float_t errSX [5];
 998     Float_t errPh [5];
 999
1000     for ( iPhi = 0; iPhi < iPhiMax; ++iPhi )
1001     {
1002       char hXtitle[80];
1003       char hYtitle[80];
1004
1005       phiMin = 72*iPhi     ;
1006       phiMax = 72*iPhi + 72;
1007
1008       TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1009
1010       if ( firstChamber == lastChamber ) {
1011         sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d",phiMin,phiMax,firstChamber);
1012         sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d",phiMin,phiMax,lastChamber);
1013       }
1014       else{
1015         sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
1016         sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
1017       }
1018
1019       hDeltaY [iPhi] -> SetTitle(hYtitle);
1020       hDeltaX [iPhi] -> SetTitle(hXtitle);
1021
1022       hDeltaY [iPhi]      -> Fit("fY","R","E");
1023       sigmaY  [iPhi] = fY -> GetParameter(2)  ;
1024       errSY   [iPhi] = fY -> GetParError(2)   ;
1025
1026       rmsX    [iPhi] = hDeltaX [iPhi] -> GetRMS();
1027       errSX   [iPhi] = hDeltaX [iPhi] -> GetRMSError();
1028
1029       phi     [iPhi] = 72*iPhi + 36 ;
1030       errPh   [iPhi] = 36;
1031     }
1032
1033     //--------------------------------------------------------------------------------------------------------------
1034     //For plotting resolution in function of the angle of incidence
1035
1036     TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1037     c1-> Divide(1,2);
1038     c1->cd(1);
1039
1040     TGraphErrors * GraphX = new TGraphErrors( iPhiMax, phi, rmsX, errPh, errSX);
1041     GraphX->SetTitle("Resolution en X (Phi)");
1042     GraphX->Draw("ALP");
1043
1044     c1->cd(2);
1045     TGraphErrors * GraphY = new TGraphErrors( iPhiMax, phi, sigmaY, errPh, errSY);
1046     GraphY->SetTitle("Resolution en Y (Phi)");
1047     GraphY->Draw("ALP");
1048
1049     cout<<"\n\n\n";
1050
1051     MUONLoader->UnloadTracks();
1052
1053 }
1054
1055
1056
1057
1058
1059
1060
1061 /*****************************************************************************************************************/
1062 /*****************************************************************************************************************/
1063                                        /*RESOLUTION IN FUNCTION OF THETA*/
1064
1065 void resolutionTheta( Int_t event2Check=0, char * filename="galice.root" )
1066 {
1067   cout<<"\nChamber(s) chosen;\nFirst chamber:";
1068   cin>>firstChamber;
1069   cout<<"Last chamber:";
1070   cin>>lastChamber;
1071   cout<<"\n\n";
1072
1073   TH1F * hDeltaX[9];
1074   TH1F * hDeltaY[9];
1075
1076   hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
1077   hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
1078   hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
1079   hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
1080   hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
1081   hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
1082   hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
1083   hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
1084   hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
1085
1086   hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
1087   hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
1088   hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
1089   hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
1090   hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
1091   hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
1092   hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
1093   hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
1094   hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
1095
1096
1097 //Creating Run Loader and openning file containing Hits
1098 //--------------------------------------------------------------------------------------------------------------
1099   AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
1100   if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
1101
1102
1103 //Getting MUONLoader
1104 //--------------------------------------------------------------------------------------------------------------
1105   AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1106   MUONLoader -> LoadTracks("READ");
1107   MUONLoader -> LoadRecPoints("READ");
1108
1109
1110 //Creating a MUON data container
1111 //--------------------------------------------------------------------------------------------------------------
1112   AliMUONData muondata(MUONLoader,"MUON","MUON");
1113
1114
1115   nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1116
1117     //Loop on events
1118     for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1119     {
1120       if ( event2Check!=0 )
1121           iEvent=event2Check;
1122       printf("\r>>> Event %d ",iEvent);
1123       RunLoader->GetEvent(iEvent);
1124
1125       //Addressing
1126       muondata.SetTreeAddress("RT");
1127       muondata.GetRecTracks();
1128       tracks = muondata.RecTracks();
1129
1130       //Loop on track
1131       nTracks = (Int_t) tracks -> GetEntriesFast();
1132       for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1133       {
1134         track       = (AliMUONTrack*) tracks ->At(iTrack)          ;
1135         trackParams =                 track  ->GetTrackParamAtCluster();
1136         hits        =                 track  ->GetHitForRecAtHit() ;
1137
1138         //Loop on hits
1139         nHits = (Int_t) hits -> GetEntriesFast();
1140         for ( iHit = 0; iHit < nHits; ++iHit )
1141         {
1142           trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
1143           hit        = (AliMUONHitForRec* ) hits        -> At(iHit);
1144           chamberNbr = hit -> GetChamberNumber();
1145           if ( chamberNbr >= firstChamber -1 ) {
1146             if ( chamberNbr < lastChamber ) {
1147               //Positions
1148               Double_t traX, traY;
1149               Double_t hitX, hitY, hitZ;
1150               Double_t aveY, aveX;
1151
1152               //Angles
1153               Double_t theta;
1154               Int_t    iTheta;
1155
1156               //Resolution
1157               Double_t deltaX;
1158               Double_t deltaY;
1159
1160               traX = trackParam -> GetNonBendingCoor();
1161               traY = trackParam -> GetBendingCoor()   ;
1162
1163               hitX = hit        -> GetNonBendingCoor();
1164               hitY = hit        -> GetBendingCoor()   ;
1165               hitZ = hit        -> GetZ();
1166
1167               aveX = 1./2.*(traX + hitX);
1168               aveY = 1./2.*(traY + hitY);
1169
1170               //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
1171               theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ ));
1172
1173               iTheta = int( theta );
1174               if( theta > 10 ) iTheta = 9;
1175               if( theta < 1  ) iTheta = 1;
1176
1177               deltaX = traX - hitX;
1178               deltaY = traY - hitY;
1179
1180               hDeltaX [iTheta-1] -> Fill( deltaX );
1181               hDeltaY [iTheta-1] -> Fill( deltaY );
1182             }
1183           }
1184
1185         }
1186         //End loop on hits
1187
1188       }
1189       //End loop on tracks
1190       muondata.ResetRecTracks();
1191       if ( event2Check!=0 )
1192           iEvent=nEvents;
1193
1194     }
1195     //End loop on events
1196
1197
1198     Int_t iTheta;
1199     Int_t iThetaMax = 9;
1200     Int_t thetaMin, thetaMax;
1201
1202     Float_t theta [9];
1203     Float_t sigmaY[9];
1204     Float_t errSY [9];
1205     Float_t rmsX  [9];
1206     Float_t errSX [9];
1207     Float_t ErrTh [9];
1208
1209     for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1210     {
1211       char hXtitle[80];
1212       char hYtitle[80];
1213
1214       //To find the polar angle
1215       thetaMin = 178 - iTheta;
1216       thetaMax = 179 - iTheta;
1217
1218       TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1219
1220       if ( firstChamber == lastChamber ) {
1221         sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
1222         sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
1223       }
1224       else{
1225         sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1226         sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1227       }
1228
1229       hDeltaY [iTheta] -> SetTitle(hYtitle);
1230       hDeltaX [iTheta] -> SetTitle(hXtitle);
1231
1232       hDeltaY [iTheta]      -> Fit("fY","R","E");
1233       sigmaY  [iTheta] = fY -> GetParameter(2)  ;
1234       errSY   [iTheta] = fY -> GetParError(2)   ;
1235
1236       rmsX    [iTheta] = hDeltaX [iTheta] -> GetRMS();
1237       errSX   [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1238
1239       theta   [iTheta] = 178.5 - iTheta;
1240       ErrTh   [iTheta] = 0.5;
1241     }
1242
1243     //--------------------------------------------------------------------------------------------------------------
1244     //For plotting resolution in function of the angle of incidence
1245
1246     TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1247     c1-> Divide(1,2);
1248     c1->cd(1);
1249
1250     TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, ErrTh, errSX);
1251     GraphX->SetTitle("Resolution en X (Theta)");
1252     GraphX->Draw("ALP");
1253
1254     c1->cd(2);
1255     TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, ErrTh, errSY);
1256     GraphY->SetTitle("Resolution en Y (Theta)");
1257     GraphY->Draw("ALP");
1258
1259     cout<<"\n\n\n";
1260     MUONLoader->UnloadTracks();
1261
1262 }
1263
1264
1265
1266
1267
1268 /*****************************************************************************************************************/
1269 /*****************************************************************************************************************/
1270                               /*RESOLUTION IN FUNCTION OF THETA OF INCIDENCE*/
1271
1272 void resolutionThetaI( Int_t event2Check=0, char * filename="galice.root" )
1273 {
1274   cout<<"\nChamber(s) chosen;\nFirst chamber:";
1275   cin>>firstChamber;
1276   cout<<"Last chamber:";
1277   cin>>lastChamber;
1278   cout<<"\n\n";
1279
1280   TH1F * hDeltaX[11];
1281   TH1F * hDeltaY[11];
1282
1283   hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
1284   hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
1285   hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
1286   hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
1287   hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
1288   hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
1289   hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
1290   hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
1291   hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
1292   hDeltaX[9] = new TH1F ("hDeltaX9" , "No Interest", 100, -0.3, 0.3);
1293   hDeltaX[10]= new TH1F ("hDeltaX10", "No Interest", 100, -0.3, 0.3);
1294
1295   hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
1296   hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
1297   hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
1298   hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
1299   hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
1300   hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
1301   hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
1302   hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
1303   hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
1304   hDeltaY[9] = new TH1F ("hDeltaY9" , "No Interest", 500, -0.1, 0.1);
1305   hDeltaY[10]= new TH1F ("hDeltaY10", "No Interest", 500, -0.1, 0.1);
1306
1307
1308 //Creating Run Loader and openning file containing Hits
1309 //--------------------------------------------------------------------------------------------------------------
1310   AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
1311   if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
1312
1313
1314 //Getting MUONLoader
1315 //--------------------------------------------------------------------------------------------------------------
1316   AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1317   MUONLoader -> LoadTracks("READ");
1318   MUONLoader -> LoadRecPoints("READ");
1319
1320
1321 //Creating a MUON data container
1322 //--------------------------------------------------------------------------------------------------------------
1323   AliMUONData muondata(MUONLoader,"MUON","MUON");
1324
1325
1326   nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1327
1328     //Loop on events
1329     for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1330     {
1331       if ( event2Check!=0 )
1332           iEvent=event2Check;
1333       printf("\r>>> Event %d ",iEvent);
1334       RunLoader->GetEvent(iEvent);
1335
1336       //Addressing
1337       muondata.SetTreeAddress("RT");
1338       muondata.GetRecTracks();
1339       tracks = muondata.RecTracks();
1340
1341       //Loop on track
1342       nTracks = (Int_t) tracks -> GetEntriesFast();
1343       for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1344       {
1345         track       = (AliMUONTrack*) tracks ->At(iTrack)          ;
1346         trackParams =                 track  ->GetTrackParamAtCluster();
1347         hits        =                 track  ->GetHitForRecAtHit() ;
1348
1349         //Loop on hits
1350         nHits = (Int_t) hits -> GetEntriesFast();
1351         for ( iHit = 0; iHit < nHits; ++iHit )
1352         {
1353           trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
1354           hit        = (AliMUONHitForRec* ) hits        -> At(iHit);
1355           chamberNbr = hit -> GetChamberNumber();
1356           if ( chamberNbr >= firstChamber - 1 ) {
1357             if ( chamberNbr < lastChamber ) {
1358               //Momentum
1359               Double_t Px, Py, Pz, Pr;
1360
1361               //Positions
1362               Double_t traX, traY;
1363               Double_t hitX, hitY;
1364
1365               //Resolution
1366               Double_t deltaX;
1367               Double_t deltaY;
1368
1369               //Angle
1370               Double_t theta;
1371               Int_t    iTheta;
1372
1373               Px = trackParam -> Px()   ;
1374               Py = trackParam -> Py()   ;
1375               Pz = trackParam -> Pz()   ;
1376               Pr = sqrt( Px*Px + Py*Py );
1377
1378               traX = trackParam -> GetNonBendingCoor();
1379               traY = trackParam -> GetBendingCoor();
1380
1381               hitX = hit        -> GetNonBendingCoor();
1382               hitY = hit        -> GetBendingCoor();
1383
1384               //The calculation of theta, the angle of incidence:
1385               theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr))
1386 ;
1387               if ( theta < 79 ) iTheta = 0;
1388               else iTheta = int( theta - 79.);
1389
1390               deltaX = traX - hitX;
1391               deltaY = traY - hitY;
1392
1393               hDeltaX [iTheta] -> Fill (deltaX);
1394               hDeltaY [iTheta] -> Fill (deltaY);
1395             }
1396           }
1397
1398         }
1399         //End loop on hits
1400
1401       }
1402       //End loop on tracks
1403       muondata.ResetRecTracks();
1404       if ( event2Check!=0 )
1405           iEvent=nEvents;
1406
1407     }
1408     //End loop on events
1409
1410
1411     Int_t iTheta;
1412     Int_t iThetaMax = 11;
1413     Int_t thetaMin, thetaMax;
1414
1415     Float_t theta [11];
1416     Float_t sigmaY[11];
1417     Float_t errSY [11];
1418     Float_t rmsX  [11];
1419     Float_t errSX [11];
1420     Float_t errTh [11];
1421
1422     for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1423     {
1424       char hXtitle[80];
1425       char hYtitle[80];
1426
1427       thetaMin = iTheta + 79;
1428       thetaMax = iTheta + 80;
1429
1430       TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1431
1432       if ( firstChamber == lastChamber ) {
1433         sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
1434         sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
1435       }
1436       else{
1437         sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1438         sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1439       }
1440
1441       hDeltaY [iTheta] -> SetTitle(hYtitle);
1442       hDeltaX [iTheta] -> SetTitle(hXtitle);
1443
1444       hDeltaY [iTheta]      -> Fit("fY","R","E");
1445       sigmaY  [iTheta] = fY -> GetParameter(2)  ;
1446       errSY   [iTheta] = fY -> GetParError(2)   ;
1447
1448       rmsX    [iTheta] = hDeltaX [iTheta] -> GetRMS();
1449       errSX   [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1450
1451       theta   [iTheta] = iTheta + 79.5 ;
1452       errTh   [iTheta] = 0.5;
1453     }
1454
1455     //--------------------------------------------------------------------------------------------------------------
1456     //For plotting resolution in function of the angle of incidence
1457
1458     TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1459     c1-> Divide(1,2);
1460     c1->cd(1);
1461
1462     TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, errTh, errSX);
1463     GraphX->SetTitle("Resolution en X (Theta)");
1464     GraphX->Draw("ALP");
1465
1466     c1->cd(2);
1467     TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, errTh, errSY);
1468     GraphY->SetTitle("Resolution en Y (Theta)");
1469     GraphY->Draw("ALP");
1470
1471     cout<<"\n\n\n";
1472     MUONLoader->UnloadTracks();
1473
1474 }
1475
1476
1477
1478