MUON/MUONResoEffChamber.C

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