MUON/MUONResoEffChamber.C

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