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    3 pt cuts. Results are compared to the reference (red curves).
  21    To be used with (at least) 10000 events as follows
  22    AliGenBox * gener = new AliGenBox(1);
  23    gener->SetPtRange(0.,10.);
  24    gener->SetPhiRange(0., 360.);
  25    gener->SetThetaRange(171.000,178.001);
  26    gener->SetPart(13);           // or -13
  27    gener->SetOrigin(0.,0., 0.);
  28    gener->SetSigma(0.0, 0.0, 0.0);
  29
  30    Author: Fabien Guerin, LPC Clermont-Ferrand, Jan. 2006
  31 */
  32
  33
  34 // ROOT includes
  35 #include "TBranch.h"
  36 #include "TClonesArray.h"
  37 #include "TFile.h"
  38 #include "TF1.h"
  39 #include "TH1.h"
  40 #include "TH2.h"
  41 #include "TParticle.h"
  42 #include "TTree.h"
  43 #include "TMath.h"
  44 #include "TStyle.h"
  45 #include "TCanvas.h"
  46 #include "TLegend.h"
  47
  48
  49 // STEER includes
  50 #include "AliRun.h"
  51 #include "AliDetector.h"
  52 #include "AliRunLoader.h"
  53 #include "AliHeader.h"
  54 #include "AliLoader.h"
  55 #include "AliStack.h"
  56 #include "AliSegmentation.h"
  57 #include "AliMC.h"
  58
  59 // MUON includes
  60 #include "AliMUON.h"
  61 #include "AliMUONData.h"
  62 #include "AliMUONHit.h"
  63 #include "AliMUONChamber.h"
  64 #include "AliMUONConstants.h"
  65 #include "AliMUONDigit.h"
  66 #include "AliMUONRawCluster.h"
  67 #include "AliMUONGlobalTrigger.h"
  68 #include "AliMUONLocalTrigger.h"
  69 #include "AliMUONTrack.h"
  70
  71 Double_t fitArc(Double_t *x,Double_t *par)
  72 {
  73   Double_t h1=par[1]*(x[0]+par[2]);
  74   return par[0]*TMath::ATan(TMath::Power(h1,6))+par[3];
  75 }
  76
  77 Double_t fitArch(Double_t *x,Double_t *par)
  78 {
  79   Double_t h0=2*par[0]/(TMath::Pi());
  80   Double_t h1=par[1]*x[0]+par[2]*x[0]*x[0];
  81   return h0*TMath::TanH(h1)+par[3];
  82 }
  83
  84 void MUONTriggerEfficiencyPt(char filename[10]="galice.root")
  85 {
  86
  87 // define style
  88     TStyle *st1 = new TStyle("st1","My STYLE");
  89
  90     st1->SetOptStat(0);
  91     st1->SetOptFit(111);
  92     st1->SetFrameFillColor(10);
  93     st1->SetCanvasColor(10);
  94     st1->SetFillColor(10);
  95     st1->SetTitleBorderSize(0);
  96     st1->SetTitleOffset(1.2,"XY");
  97     st1->SetTitleSize(0.05,"XY");
  98     st1->SetLabelSize(0.045,"XY");
  99     st1->SetLabelFont(22,"XY");
 100     st1->SetTitleFont(22,"XY");
 101     st1->SetOptTitle(0);
 102     st1->SetStatColor(10);
 103     st1->SetStatFont(62);
 104     st1->SetFitFormat("4.4g");
 105     st1->SetPadLeftMargin(0.15);
 106     st1->SetPadRightMargin(0.06);
 107     st1->SetPadTopMargin(0.06);
 108     st1->SetPadBottomMargin(0.15);
 109     st1->cd();
 110
 111     gROOT->ForceStyle();
 112     //TGaxis::SetMaxDigits(3);
 113
 114 // beginning of macro
 115     TH1F *ptapt    = new TH1F("ptapt","",50,0.15,5.);
 116     TH1F *ptlpt    = new TH1F("ptlpt","",50,0.15,5.);
 117     TH1F *pthpt    = new TH1F("pthpt","",50,0.15,5.);
 118     TH1F *ptcoinc  = new TH1F("ptcoinc","",50,0.15,5.);
 119
 120     TParticle *particle;
 121     AliStack* stack;
 122
 123     Double_t coincmuon=0;
 124     Double_t aptmuon=0;
 125     Double_t lptmuon=0;
 126     Double_t hptmuon=0;
 127     Double_t ptmu;
 128
 129 // output file
 130     char digitdat[100];
 131     sprintf(digitdat,"MUONTriggerEfficiencyPt.out");
 132     FILE *fdat=fopen(digitdat,"w");
 133
 134 // Creating Run Loader and openning file containing Hits
 135     AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
 136
 137     if (RunLoader ==0x0) {
 138         printf(">>> Error : Error Opening %s file \n",filename);
 139         return;
 140     }
 141
 142     Int_t nevents = RunLoader->GetNumberOfEvents();
 143
 144 // loading hits
 145     AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader");
 146     MUONLoader->LoadHits("READ");
 147     AliMUONData data_hits(MUONLoader,"MUON","MUON");
 148
 149     TClonesArray * globalTrigger;
 150     AliMUONGlobalTrigger * gloTrg;
 151
 152 // loading trigger
 153     MUONLoader->LoadDigits("READ");
 154     AliMUONData muondata(MUONLoader,"MUON","MUON");
 155
 156 // Loading kine
 157     RunLoader->LoadKinematics("READ");
 158
 159     for (Int_t ievent=0; ievent<nevents; ievent++) {  // Event loop
 160
 161         RunLoader->GetEvent(ievent);
 162
 163         if (ievent%500==0) printf("ievent = %d \n",ievent);
 164
 165
 166 // kine
 167         Int_t iparticle, nparticles;
 168         stack = RunLoader->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         muondata.SetTreeAddress("D,GLT");
 179         muondata.GetTriggerD();
 180
 181         globalTrigger = muondata.GlobalTrigger();
 182
 183         Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1
 184         for (Int_t iglobal=0; iglobal<nglobals; iglobal++) { // Global Trigger
 185             gloTrg = static_cast<AliMUONGlobalTrigger*>(globalTrigger->At(iglobal));
 186
 187             if (gloTrg->SinglePlusApt()>=1 || gloTrg->SingleMinusApt()>=1 || gloTrg->SingleUndefApt()) {
 188                 aptmuon++;
 189                 ptapt->Fill(ptmu);
 190             }
 191             if (gloTrg->SinglePlusLpt()>=1 || gloTrg->SingleMinusLpt()>=1 || gloTrg->SingleUndefLpt()) {
 192                 lptmuon++;
 193                 ptlpt->Fill(ptmu);
 194             }
 195             if (gloTrg->SinglePlusHpt()>=1 || gloTrg->SingleMinusHpt()>=1 || gloTrg->SingleUndefHpt()) {
 196                 hptmuon++;
 197                 pthpt->Fill(ptmu);
 198             }
 199         } // end of loop on Global Trigger
 200         muondata.ResetTrigger();
 201
 202 // Hits
 203         RunLoader->GetEvent(ievent);
 204         data_hits.SetTreeAddress("H");
 205
 206         Int_t itrack, ntracks, NbHits[4];
 207         Int_t SumNbHits;
 208
 209         for (Int_t j=0; j<4; j++) NbHits[j]=0;
 210
 211         ntracks = (Int_t) data_hits.GetNtracks();
 212
 213         for (itrack=0; itrack<ntracks; itrack++) { // track loop
 214          data_hits.GetTrack(itrack);
 215
 216          Int_t ihit, nhits;
 217          nhits = (Int_t) data_hits.Hits()->GetEntriesFast();
 218          AliMUONHit* mHit;
 219
 220           for (ihit=0; ihit<nhits; ihit++) {
 221             mHit = static_cast<AliMUONHit*>(data_hits.Hits()->At(ihit));
 222             Int_t Nch        = mHit->Chamber();
 223             Int_t hittrack   = mHit->Track();
 224             Int_t IdPart     = mHit->Particle();
 225
 226             for (Int_t j=0;j<4;j++) {
 227               Int_t kch=11+j;
 228               if (hittrack==0) {
 229                 if (Nch==kch && (IdPart==-13 || IdPart==13) && NbHits[j]==0) NbHits[j]++;
 230               }
 231             }
 232           }
 233           data_hits.ResetHits();
 234         } // end track loop
 235
 236
 237 // 3/4 coincidence
 238         SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3];
 239
 240         if (SumNbHits==3 || SumNbHits==4) {
 241             coincmuon++;
 242             ptcoinc->Fill(ptmu);
 243         }
 244
 245       } // end loop on event
 246
 247       MUONLoader->UnloadHits();
 248       MUONLoader->UnloadDigits();
 249       MUONLoader->UnloadRecPoints();
 250       RunLoader->UnloadKinematics();
 251
 252       delete RunLoader;
 253
 254       fprintf(fdat,"\n");
 255       fprintf(fdat,"\n");
 256       fprintf(fdat," Number of events = %d \n",nevents);
 257       fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon);
 258       fprintf(fdat," Number of dimuons with 3/4 coinc Apt cut = %d \n",(Int_t)aptmuon);
 259       fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon);
 260       fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon);
 261       fprintf(fdat,"\n");
 262
 263       Double_t efficiency,error;
 264
 265       efficiency=aptmuon/coincmuon;
 266       error=efficiency*TMath::Sqrt((aptmuon+coincmuon)/(aptmuon*coincmuon));
 267       fprintf(fdat," Efficiency Apt cut = %4.4f +/- %4.4f\n",efficiency,error);
 268
 269       efficiency=lptmuon/coincmuon;
 270       error=efficiency*TMath::Sqrt((lptmuon+coincmuon)/(lptmuon*coincmuon));
 271       fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error);
 272
 273       efficiency=hptmuon/coincmuon;
 274       error=efficiency*TMath::Sqrt((hptmuon+coincmuon)/(hptmuon*coincmuon));
 275       fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error);
 276
 277       fclose(fdat);
 278
 279       Double_t x1,x2,xval,xerr;
 280
 281       for (Int_t i=0;i<50;i++) {
 282           x1=ptcoinc->GetBinContent(i+1);
 283
 284           x2=ptapt->GetBinContent(i+1);
 285           if (x1!=0 && x2!=0) {
 286               xval=x2/x1;
 287               xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);
 288               ptapt->SetBinContent(i+1,xval);
 289               ptapt->SetBinError(i+1,0);
 290           } else {
 291               ptapt->SetBinContent(i+1,0.);
 292               ptapt->SetBinError(i+1,0.);
 293           }
 294
 295           x2=ptlpt->GetBinContent(i+1);
 296           if (x1!=0 && x2!=0) {
 297               xval=x2/x1;
 298               xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);
 299               ptlpt->SetBinContent(i+1,xval);
 300               ptlpt->SetBinError(i+1,0);
 301           } else {
 302               ptlpt->SetBinContent(i+1,0.);
 303               ptlpt->SetBinError(i+1,0.);
 304           }
 305
 306           x2=pthpt->GetBinContent(i+1);
 307           if (x1!=0 && x2!=0) {
 308               xval=x2/x1;
 309               xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);
 310               pthpt->SetBinContent(i+1,xval);
 311               pthpt->SetBinError(i+1,0);
 312           } else {
 313               pthpt->SetBinContent(i+1,0.);
 314               pthpt->SetBinError(i+1,0.);
 315           }
 316       }
 317
 318
 319       TF1 *fitapt = new TF1("fitapt",fitArch,0.,5.,4);
 320       TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4);
 321       TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4);
 322
 323       TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400);
 324       c1->Divide(3,1);
 325
 326       c1->cd(1);
 327       ptapt->SetMinimum(0.);
 328       ptapt->SetMaximum(1.05);
 329       ptapt->SetTitle("");
 330       ptapt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 331       ptapt->GetYaxis()->SetTitle("Efficiency");
 332       //ptapt->GetXaxis()->SetRange(3);
 333       ptapt->Draw("");
 334       fitapt->SetParameters(1.853,1.57,-0.0203,-0.178);
 335       fitapt->SetLineColor(2);
 336       fitapt->SetLineWidth(2);
 337       fitapt->Draw("SAME");
 338       TLegend * leg = new TLegend(0.5,0.38,0.85,0.53);
 339       leg->SetBorderSize(0);
 340       leg->SetFillColor(0);
 341       leg->SetTextSize(0.05);
 342       leg->SetTextFont(22);
 343       leg->SetHeader("Apt trigger pt cut");
 344       leg->AddEntry(fitapt,"Ref.","l");
 345       leg->AddEntry(ptapt,"New","l");
 346       leg->Draw("SAME");
 347
 348       c1->cd(2);
 349       ptlpt->SetMinimum(0.);
 350       ptlpt->SetMaximum(1.05);
 351       ptlpt->SetTitle("");
 352       ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 353       ptlpt->GetYaxis()->SetTitle("Efficiency");
 354       //ptlpt->GetXaxis()->SetRange(3);
 355       ptlpt->Draw("");
 356       fitlpt->SetParameters(0.602,0.774,0.273,0.048);
 357       fitlpt->SetLineColor(2);
 358       fitlpt->SetLineWidth(2);
 359       fitlpt->Draw("SAME");
 360       leg = new TLegend(0.5,0.38,0.85,0.53);
 361       leg->SetBorderSize(0);
 362       leg->SetFillColor(0);
 363       leg->SetTextSize(0.05);
 364       leg->SetTextFont(22);
 365       leg->SetHeader("Lpt trigger pt cut");
 366       leg->AddEntry(fitlpt,"Ref.","l");
 367       leg->AddEntry(ptlpt,"New","l");
 368       leg->Draw("SAME");
 369
 370       c1->cd(3);
 371       pthpt->SetMinimum(0.);
 372       pthpt->SetMaximum(1.05);
 373       pthpt->SetTitle("");
 374       pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 375       pthpt->GetYaxis()->SetTitle("Efficiency");
 376       //pthpt->GetXaxis()->SetRange(3);
 377       pthpt->Draw("");
 378       fithpt->SetParameters(0.627,0.393,0.855,0.0081);
 379       fithpt->SetLineColor(2);
 380       fithpt->SetLineWidth(2);
 381       fithpt->Draw("SAME");
 382       leg = new TLegend(0.5,0.38,0.85,0.53);
 383       leg->SetBorderSize(0);
 384       leg->SetFillColor(0);
 385       leg->SetTextSize(0.05);
 386       leg->SetTextFont(22);
 387       leg->SetHeader("Hpt trigger pt cut");
 388       leg->AddEntry(fithpt,"Ref.","l");
 389       leg->AddEntry(pthpt,"New","l");
 390       leg->Draw("SAME");
 391
 392       c1->SaveAs("MUONTriggerEfficiencyPt.gif");
 393       c1->SaveAs("MUONTriggerEfficiencyPt.eps");
 394
 395 }