MUON/MUONTriggerEfficiencyPt.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
  19 // Macro to produce trigger single muon efficiency versus pt plots for the
  20 // 2 pt cuts.
  21 // see full description in the README file
  22 // Author: Fabien Guerin (LPC)
  23
  24 // ROOT includes
  25 #include "TBranch.h"
  26 #include "TClonesArray.h"
  27 #include "TFile.h"
  28 #include "TF1.h"
  29 #include "TH1.h"
  30 #include "TH2.h"
  31 #include "TParticle.h"
  32 #include "TTree.h"
  33 #include "TMath.h"
  34 #include "TStyle.h"
  35 #include "TCanvas.h"
  36 #include "TLegend.h"
  37 #include "Riostream.h"
  38
  39 // STEER includes
  40 #include "AliRun.h"
  41 #include "AliRunLoader.h"
  42 #include "AliHeader.h"
  43 #include "AliLoader.h"
  44 #include "AliStack.h"
  45
  46 // MUON includes
  47 #include "AliMUONSimData.h"
  48 #include "AliMUONRecData.h"
  49 #include "AliMUONHit.h"
  50 #include "AliMUONDigit.h"
  51 #include "AliMUONGlobalTrigger.h"
  52 #include "AliMUONLocalTrigger.h"
  53 #include "AliMUONTrack.h"
  54
  55 Double_t fitArc(Double_t *x,Double_t *par)
  56 {
  57   Double_t h1=par[1]*(x[0]+par[2]);
  58   return par[0]*TMath::ATan(TMath::Power(h1,6))+par[3];
  59 }
  60
  61 Double_t fitArch(Double_t *x,Double_t *par)
  62 {
  63   Double_t h0=2*par[0]/(TMath::Pi());
  64   Double_t h1=par[1]*x[0]+par[2]*x[0]*x[0];
  65   return h0*TMath::TanH(h1)+par[3];
  66 }
  67
  68 void MUONTriggerEfficiencyPt(char filenameSim[10]="galice.root",
  69                              char filenameRec[10]="galice.root",
  70                              Bool_t readFromRP = 0)
  71 {
  72
  73 // define style
  74     TStyle *st1 = new TStyle("st1","My STYLE");
  75
  76     st1->SetOptStat(0);
  77     st1->SetOptFit(111);
  78     st1->SetFrameFillColor(10);
  79     st1->SetCanvasColor(10);
  80     st1->SetFillColor(10);
  81     st1->SetTitleBorderSize(0);
  82     st1->SetTitleOffset(1.2,"XY");
  83     st1->SetTitleSize(0.05,"XY");
  84     st1->SetLabelSize(0.045,"XY");
  85     st1->SetLabelFont(22,"XY");
  86     st1->SetTitleFont(22,"XY");
  87     st1->SetOptTitle(0);
  88     st1->SetStatColor(10);
  89     st1->SetStatFont(62);
  90     st1->SetFitFormat("4.4g");
  91     st1->SetPadLeftMargin(0.15);
  92     st1->SetPadRightMargin(0.06);
  93     st1->SetPadTopMargin(0.06);
  94     st1->SetPadBottomMargin(0.15);
  95     st1->cd();
  96
  97 //    gROOT->ForceStyle();
  98     //TGaxis::SetMaxDigits(3);
  99
 100 // beginning of macro
 101     TH1F *ptlpt    = new TH1F("ptlpt","",50,0.15,5.);
 102     TH1F *pthpt    = new TH1F("pthpt","",50,0.15,5.);
 103     TH1F *ptcoinc  = new TH1F("ptcoinc","",50,0.15,5.);
 104
 105     TParticle *particle;
 106     AliStack* stack;
 107
 108     Int_t nevents;
 109     Double_t coincmuon=0;
 110     Double_t lptmuon=0;
 111     Double_t hptmuon=0;
 112     Double_t ptmu=0;
 113
 114 // output file
 115     char digitdat[100];
 116     char currentfile[100];
 117     sprintf(digitdat,"MUONTriggerEfficiencyPt.out");
 118     FILE *fdat=fopen(digitdat,"w");
 119
 120 // Initialise AliRoot
 121 // Creating Run Loader and openning file containing Hits
 122     AliRunLoader * RunLoaderSim = AliRunLoader::Open(filenameSim,"MUONFolderSim","READ");
 123
 124     if (RunLoaderSim ==0x0) {
 125         printf(">>> Error : Error Opening %s file \n",currentfile);
 126         return;
 127     }
 128
 129     AliRunLoader * RunLoaderRec = AliRunLoader::Open(filenameRec,"MUONFolder","READ");
 130
 131     if (RunLoaderRec ==0x0) {
 132         printf(">>> Error : Error Opening %s file \n",currentfile);
 133         return;
 134     }
 135
 136     nevents = RunLoaderSim->GetNumberOfEvents();
 137
 138     AliLoader * MUONLoaderSim = RunLoaderSim->GetLoader("MUONLoader");
 139     AliLoader * MUONLoaderRec = RunLoaderRec->GetLoader("MUONLoader");
 140
 141     if (!readFromRP) {
 142         cout << " reading from digits \n";
 143         MUONLoaderSim->LoadDigits("READ");
 144     } else {
 145         cout << " reading from RecPoints \n";
 146         MUONLoaderRec->LoadRecPoints("READ");
 147    }
 148     MUONLoaderSim->LoadHits("READ");
 149     RunLoaderSim->LoadKinematics("READ");
 150
 151
 152 // Creating MUON data containers
 153     AliMUONSimData muondataSim(MUONLoaderSim,"MUON","MUON");
 154     AliMUONRecData muondataRec(MUONLoaderRec,"MUON","MUON");
 155
 156     TClonesArray * globalTrigger;
 157     AliMUONGlobalTrigger * gloTrg;
 158
 159     for (Int_t ievent=0; ievent<nevents; ievent++) {  // Event loop
 160
 161         RunLoaderSim->GetEvent(ievent);
 162         RunLoaderRec->GetEvent(ievent);
 163
 164         if (ievent%500==0) printf("ievent = %d \n",ievent);
 165
 166 // kine
 167         Int_t iparticle, nparticles;
 168         stack = RunLoaderSim->Stack();
 169         nparticles = (Int_t) stack->GetNtrack();
 170         for (iparticle=0; iparticle<nparticles; iparticle++) {
 171             particle = stack->Particle(iparticle);
 172             Float_t pt = particle->Pt();
 173             Int_t pdgcode = particle->GetPdgCode();
 174             if (pdgcode==-13 || pdgcode==13) ptmu = pt;
 175         }
 176
 177 // trigger
 178         if (!readFromRP) {
 179             muondataSim.SetTreeAddress("D,GLT");
 180             muondataSim.GetTriggerD();
 181             globalTrigger = muondataSim.GlobalTrigger();
 182         } else {
 183             muondataRec.SetTreeAddress("RC,TC");
 184             muondataRec.GetTrigger();
 185             globalTrigger = muondataRec.GlobalTrigger();
 186
 187         }
 188
 189         Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1
 190         for (Int_t iglobal=0; iglobal<nglobals; iglobal++) { // Global Trigger
 191             gloTrg = static_cast<AliMUONGlobalTrigger*>(globalTrigger->At(iglobal));
 192
 193             if (gloTrg->SingleLpt()>=1 ) {
 194                 lptmuon++;
 195                 ptlpt->Fill(ptmu);
 196             }
 197             if (gloTrg->SingleHpt()>=1 ) {
 198                 hptmuon++;
 199                 pthpt->Fill(ptmu);
 200             }
 201         } // end of loop on Global Trigger
 202         muondataSim.ResetTrigger();
 203         muondataRec.ResetTrigger();
 204
 205 // Hits
 206         muondataSim.SetTreeAddress("H");
 207
 208         Int_t itrack, ntracks, NbHits[4];
 209         Int_t SumNbHits;
 210
 211         for (Int_t j=0; j<4; j++) NbHits[j]=0;
 212
 213         ntracks = (Int_t) muondataSim.GetNtracks();
 214
 215         for (itrack=0; itrack<ntracks; itrack++) { // track loop
 216             muondataSim.GetTrack(itrack);
 217
 218          Int_t ihit, nhits;
 219          nhits = (Int_t) muondataSim.Hits()->GetEntriesFast();
 220          AliMUONHit* mHit;
 221
 222           for (ihit=0; ihit<nhits; ihit++) {
 223               mHit = static_cast<AliMUONHit*>(muondataSim.Hits()->At(ihit));
 224             Int_t Nch        = mHit->Chamber();
 225             Int_t hittrack   = mHit->Track();
 226             Float_t IdPart   = mHit->Particle();
 227
 228             for (Int_t j=0;j<4;j++) {
 229               Int_t kch=11+j;
 230               if (hittrack==0) {
 231                 if (Nch==kch && (IdPart==-13 || IdPart==13) && NbHits[j]==0) NbHits[j]++;
 232                 }
 233             }
 234           }
 235           muondataSim.ResetHits();
 236         } // end track loop
 237
 238 // 3/4 coincidence
 239         SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3];
 240
 241         if (SumNbHits==3 || SumNbHits==4) {
 242             coincmuon++;
 243             ptcoinc->Fill(ptmu);
 244         }
 245
 246     } // end loop on event
 247
 248     if (coincmuon==0) {
 249         cout << " >>> <E> coincmuon = 0 after event loop " << "\n";
 250         cout << " >>> this probably means that input does not contain one (and only one) " << "\n";
 251         cout << " >>> muon track per event as it should " << "\n";
 252         cout << " >>> see README for further information " << "\n";
 253         cout << " >>> exiting now ! " << "\n";
 254         return;
 255     }
 256
 257
 258        MUONLoaderSim->UnloadHits();
 259        if (!readFromRP) {
 260            MUONLoaderSim->UnloadDigits();
 261        } else {
 262            MUONLoaderRec->UnloadRecPoints();
 263        }
 264        RunLoaderSim->UnloadKinematics();
 265
 266       fprintf(fdat,"\n");
 267       fprintf(fdat,"\n");
 268       fprintf(fdat," Number of events = %d \n",nevents);
 269       fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon);
 270       fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon);
 271       fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon);
 272       fprintf(fdat,"\n");
 273
 274       Double_t efficiency,error;
 275
 276       efficiency=lptmuon/coincmuon;
 277       error=efficiency*TMath::Sqrt((lptmuon+coincmuon)/(lptmuon*coincmuon));
 278       fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error);
 279
 280       efficiency=hptmuon/coincmuon;
 281       error=efficiency*TMath::Sqrt((hptmuon+coincmuon)/(hptmuon*coincmuon));
 282       fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error);
 283
 284       fclose(fdat);
 285
 286       Double_t x1,x2,xval,xerr;
 287
 288       for (Int_t i=0;i<50;i++) {
 289           x1=ptcoinc->GetBinContent(i+1);
 290
 291           x2=ptlpt->GetBinContent(i+1);
 292           if (x1!=0 && x2!=0) {
 293               xval=x2/x1;
 294               xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);
 295               ptlpt->SetBinContent(i+1,xval);
 296               ptlpt->SetBinError(i+1,0);
 297           } else {
 298               ptlpt->SetBinContent(i+1,0.);
 299               ptlpt->SetBinError(i+1,0.);
 300           }
 301
 302           x2=pthpt->GetBinContent(i+1);
 303           if (x1!=0 && x2!=0) {
 304               xval=x2/x1;
 305               xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);
 306               pthpt->SetBinContent(i+1,xval);
 307               pthpt->SetBinError(i+1,0);
 308           } else {
 309               pthpt->SetBinContent(i+1,0.);
 310               pthpt->SetBinError(i+1,0.);
 311           }
 312       }
 313
 314       TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4);
 315       TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4);
 316
 317       TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400);
 318       c1->Divide(2,1);
 319
 320       c1->cd(1);
 321       ptlpt->SetMinimum(0.);
 322       ptlpt->SetMaximum(1.05);
 323       ptlpt->SetTitle("");
 324       ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 325       ptlpt->GetYaxis()->SetTitle("Efficiency");
 326       //ptlpt->GetXaxis()->SetRange(3);
 327       ptlpt->Draw("");
 328       fitlpt->SetParameters(0.602,0.774,0.273,0.048);
 329       fitlpt->SetLineColor(2);
 330       fitlpt->SetLineWidth(2);
 331       fitlpt->Draw("SAME");
 332       TLegend * leg = new TLegend(0.5,0.38,0.85,0.53);
 333       leg = new TLegend(0.5,0.38,0.85,0.53);
 334       leg->SetBorderSize(0);
 335       leg->SetFillColor(0);
 336       leg->SetTextSize(0.05);
 337       leg->SetTextFont(22);
 338       leg->SetHeader("Lpt trigger pt cut");
 339       leg->AddEntry(fitlpt,"Ref.","l");
 340       leg->AddEntry(ptlpt,"New","l");
 341       leg->Draw("SAME");
 342
 343       c1->cd(2);
 344       pthpt->SetMinimum(0.);
 345       pthpt->SetMaximum(1.05);
 346       pthpt->SetTitle("");
 347       pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 348       pthpt->GetYaxis()->SetTitle("Efficiency");
 349       //pthpt->GetXaxis()->SetRange(3);
 350       pthpt->Draw("");
 351       fithpt->SetParameters(0.627,0.393,0.855,0.0081);
 352       fithpt->SetLineColor(2);
 353       fithpt->SetLineWidth(2);
 354       fithpt->Draw("SAME");
 355       leg = new TLegend(0.5,0.38,0.85,0.53);
 356       leg->SetBorderSize(0);
 357       leg->SetFillColor(0);
 358       leg->SetTextSize(0.05);
 359       leg->SetTextFont(22);
 360       leg->SetHeader("Hpt trigger pt cut");
 361       leg->AddEntry(fithpt,"Ref.","l");
 362       leg->AddEntry(pthpt,"New","l");
 363       leg->Draw("SAME");
 364
 365       c1->SaveAs("MUONTriggerEfficiencyPt.gif");
 366       c1->SaveAs("MUONTriggerEfficiencyPt.eps");
 367 }