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