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