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