From 5707a4638f6eb835dd085c5158bdc1c9f8de178e Mon Sep 17 00:00:00 2001 From: pcrochet Date: Sat, 21 Jan 2006 17:12:31 +0000 Subject: [PATCH] macro to get single muon trigger effiency versus pt (Fabien) --- MUON/MUONTriggerEfficiencyPt.C | 395 +++++++++++++++++++++++++++++++++ 1 file changed, 395 insertions(+) create mode 100644 MUON/MUONTriggerEfficiencyPt.C diff --git a/MUON/MUONTriggerEfficiencyPt.C b/MUON/MUONTriggerEfficiencyPt.C new file mode 100644 index 00000000000..293387297e5 --- /dev/null +++ b/MUON/MUONTriggerEfficiencyPt.C @@ -0,0 +1,395 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + + +/* Macro to produce trigger single muon efficiency versus pt plots for the + 3 pt cuts. Results are compared to the reference (red curves). + To be used with (at least) 10000 events as follows + AliGenBox * gener = new AliGenBox(1); + gener->SetPtRange(0.,10.); + gener->SetPhiRange(0., 360.); + gener->SetThetaRange(171.000,178.001); + gener->SetPart(13); // or -13 + gener->SetOrigin(0.,0., 0.); + gener->SetSigma(0.0, 0.0, 0.0); + + Author: Fabien Guerin, LPC Clermont-Ferrand, Jan. 2006 +*/ + + +// ROOT includes +#include "TBranch.h" +#include "TClonesArray.h" +#include "TFile.h" +#include "TF1.h" +#include "TH1.h" +#include "TH2.h" +#include "TParticle.h" +#include "TTree.h" +#include "TMath.h" +#include "TStyle.h" +#include "TCanvas.h" +#include "TLegend.h" + + +// STEER includes +#include "AliRun.h" +#include "AliDetector.h" +#include "AliRunLoader.h" +#include "AliHeader.h" +#include "AliLoader.h" +#include "AliStack.h" +#include "AliSegmentation.h" +#include "AliMC.h" + +// MUON includes +#include "AliMUON.h" +#include "AliMUONData.h" +#include "AliMUONHit.h" +#include "AliMUONChamber.h" +#include "AliMUONConstants.h" +#include "AliMUONDigit.h" +#include "AliMUONRawCluster.h" +#include "AliMUONGlobalTrigger.h" +#include "AliMUONLocalTrigger.h" +#include "AliMUONTrack.h" + +Double_t fitArc(Double_t *x,Double_t *par) +{ + Double_t h1=par[1]*(x[0]+par[2]); + return par[0]*TMath::ATan(TMath::Power(h1,6))+par[3]; +} + +Double_t fitArch(Double_t *x,Double_t *par) +{ + Double_t h0=2*par[0]/(TMath::Pi()); + Double_t h1=par[1]*x[0]+par[2]*x[0]*x[0]; + return h0*TMath::TanH(h1)+par[3]; +} + +void MUONTriggerEfficiencyPt(char filename[10]="galice.root") +{ + +// define style + TStyle *st1 = new TStyle("st1","My STYLE"); + + st1->SetOptStat(0); + st1->SetOptFit(111); + st1->SetFrameFillColor(10); + st1->SetCanvasColor(10); + st1->SetFillColor(10); + st1->SetTitleBorderSize(0); + st1->SetTitleOffset(1.2,"XY"); + st1->SetTitleSize(0.05,"XY"); + st1->SetLabelSize(0.045,"XY"); + st1->SetLabelFont(22,"XY"); + st1->SetTitleFont(22,"XY"); + st1->SetOptTitle(0); + st1->SetStatColor(10); + st1->SetStatFont(62); + st1->SetFitFormat("4.4g"); + st1->SetPadLeftMargin(0.15); + st1->SetPadRightMargin(0.06); + st1->SetPadTopMargin(0.06); + st1->SetPadBottomMargin(0.15); + st1->cd(); + + gROOT->ForceStyle(); + //TGaxis::SetMaxDigits(3); + +// beginning of macro + TH1F *ptapt = new TH1F("ptapt","",50,0.15,5.); + TH1F *ptlpt = new TH1F("ptlpt","",50,0.15,5.); + TH1F *pthpt = new TH1F("pthpt","",50,0.15,5.); + TH1F *ptcoinc = new TH1F("ptcoinc","",50,0.15,5.); + + TParticle *particle; + AliStack* stack; + + Double_t coincmuon=0; + Double_t aptmuon=0; + Double_t lptmuon=0; + Double_t hptmuon=0; + Double_t ptmu; + +// output file + char digitdat[100]; + sprintf(digitdat,"MUONTriggerEfficiencyPt.out"); + FILE *fdat=fopen(digitdat,"w"); + +// Creating Run Loader and openning file containing Hits + AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ"); + + if (RunLoader ==0x0) { + printf(">>> Error : Error Opening %s file \n",filename); + return; + } + + Int_t nevents = RunLoader->GetNumberOfEvents(); + +// loading hits + AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader"); + MUONLoader->LoadHits("READ"); + AliMUONData data_hits(MUONLoader,"MUON","MUON"); + + TClonesArray * globalTrigger; + AliMUONGlobalTrigger * gloTrg; + +// loading trigger + MUONLoader->LoadDigits("READ"); + AliMUONData muondata(MUONLoader,"MUON","MUON"); + +// Loading kine + RunLoader->LoadKinematics("READ"); + + for (Int_t ievent=0; ieventGetEvent(ievent); + + if (ievent%500==0) printf("ievent = %d \n",ievent); + + +// kine + Int_t iparticle, nparticles; + stack = RunLoader->Stack(); + nparticles = (Int_t) stack->GetNtrack(); + for (iparticle=0; iparticleParticle(iparticle); + Float_t pt = particle->Pt(); + Int_t pdgcode = particle->GetPdgCode(); + if (pdgcode==-13 || pdgcode==13) ptmu = pt; + } + +// trigger + muondata.SetTreeAddress("D,GLT"); + muondata.GetTriggerD(); + + globalTrigger = muondata.GlobalTrigger(); + + Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1 + for (Int_t iglobal=0; iglobal(globalTrigger->At(iglobal)); + + if (gloTrg->SinglePlusApt()>=1 || gloTrg->SingleMinusApt()>=1 || gloTrg->SingleUndefApt()) { + aptmuon++; + ptapt->Fill(ptmu); + } + if (gloTrg->SinglePlusLpt()>=1 || gloTrg->SingleMinusLpt()>=1 || gloTrg->SingleUndefLpt()) { + lptmuon++; + ptlpt->Fill(ptmu); + } + if (gloTrg->SinglePlusHpt()>=1 || gloTrg->SingleMinusHpt()>=1 || gloTrg->SingleUndefHpt()) { + hptmuon++; + pthpt->Fill(ptmu); + } + } // end of loop on Global Trigger + muondata.ResetTrigger(); + +// Hits + RunLoader->GetEvent(ievent); + data_hits.SetTreeAddress("H"); + + Int_t itrack, ntracks, NbHits[4]; + Int_t SumNbHits; + + for (Int_t j=0; j<4; j++) NbHits[j]=0; + + ntracks = (Int_t) data_hits.GetNtracks(); + + for (itrack=0; itrackGetEntriesFast(); + AliMUONHit* mHit; + + for (ihit=0; ihit(data_hits.Hits()->At(ihit)); + Int_t Nch = mHit->Chamber(); + Int_t hittrack = mHit->Track(); + Int_t IdPart = mHit->Particle(); + + for (Int_t j=0;j<4;j++) { + Int_t kch=11+j; + if (hittrack==0) { + if (Nch==kch && (IdPart==-13 || IdPart==13) && NbHits[j]==0) NbHits[j]++; + } + } + } + data_hits.ResetHits(); + } // end track loop + + +// 3/4 coincidence + SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3]; + + if (SumNbHits==3 || SumNbHits==4) { + coincmuon++; + ptcoinc->Fill(ptmu); + } + + } // end loop on event + + MUONLoader->UnloadHits(); + MUONLoader->UnloadDigits(); + MUONLoader->UnloadRecPoints(); + RunLoader->UnloadKinematics(); + + delete RunLoader; + + fprintf(fdat,"\n"); + fprintf(fdat,"\n"); + fprintf(fdat," Number of events = %d \n",nevents); + fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon); + fprintf(fdat," Number of dimuons with 3/4 coinc Apt cut = %d \n",(Int_t)aptmuon); + fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon); + fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon); + fprintf(fdat,"\n"); + + Double_t efficiency,error; + + efficiency=aptmuon/coincmuon; + error=efficiency*TMath::Sqrt((aptmuon+coincmuon)/(aptmuon*coincmuon)); + fprintf(fdat," Efficiency Apt cut = %4.4f +/- %4.4f\n",efficiency,error); + + efficiency=lptmuon/coincmuon; + error=efficiency*TMath::Sqrt((lptmuon+coincmuon)/(lptmuon*coincmuon)); + fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error); + + efficiency=hptmuon/coincmuon; + error=efficiency*TMath::Sqrt((hptmuon+coincmuon)/(hptmuon*coincmuon)); + fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error); + + fclose(fdat); + + Double_t x1,x2,xval,xerr; + + for (Int_t i=0;i<50;i++) { + x1=ptcoinc->GetBinContent(i+1); + + x2=ptapt->GetBinContent(i+1); + if (x1!=0 && x2!=0) { + xval=x2/x1; + xerr=xval*TMath::Sqrt((x1+x2)/x1*x2); + ptapt->SetBinContent(i+1,xval); + ptapt->SetBinError(i+1,0); + } else { + ptapt->SetBinContent(i+1,0.); + ptapt->SetBinError(i+1,0.); + } + + x2=ptlpt->GetBinContent(i+1); + if (x1!=0 && x2!=0) { + xval=x2/x1; + xerr=xval*TMath::Sqrt((x1+x2)/x1*x2); + ptlpt->SetBinContent(i+1,xval); + ptlpt->SetBinError(i+1,0); + } else { + ptlpt->SetBinContent(i+1,0.); + ptlpt->SetBinError(i+1,0.); + } + + x2=pthpt->GetBinContent(i+1); + if (x1!=0 && x2!=0) { + xval=x2/x1; + xerr=xval*TMath::Sqrt((x1+x2)/x1*x2); + pthpt->SetBinContent(i+1,xval); + pthpt->SetBinError(i+1,0); + } else { + pthpt->SetBinContent(i+1,0.); + pthpt->SetBinError(i+1,0.); + } + } + + + TF1 *fitapt = new TF1("fitapt",fitArch,0.,5.,4); + TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4); + TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4); + + TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400); + c1->Divide(3,1); + + c1->cd(1); + ptapt->SetMinimum(0.); + ptapt->SetMaximum(1.05); + ptapt->SetTitle(""); + ptapt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); + ptapt->GetYaxis()->SetTitle("Efficiency"); + //ptapt->GetXaxis()->SetRange(3); + ptapt->Draw(""); + fitapt->SetParameters(1.853,1.57,-0.0203,-0.178); + fitapt->SetLineColor(2); + fitapt->SetLineWidth(2); + fitapt->Draw("SAME"); + TLegend * leg = new TLegend(0.5,0.38,0.85,0.53); + leg->SetBorderSize(0); + leg->SetFillColor(0); + leg->SetTextSize(0.05); + leg->SetTextFont(22); + leg->SetHeader("Apt trigger pt cut"); + leg->AddEntry(fitapt,"Ref.","l"); + leg->AddEntry(ptapt,"New","l"); + leg->Draw("SAME"); + + c1->cd(2); + ptlpt->SetMinimum(0.); + ptlpt->SetMaximum(1.05); + ptlpt->SetTitle(""); + ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); + ptlpt->GetYaxis()->SetTitle("Efficiency"); + //ptlpt->GetXaxis()->SetRange(3); + ptlpt->Draw(""); + fitlpt->SetParameters(0.602,0.774,0.273,0.048); + fitlpt->SetLineColor(2); + fitlpt->SetLineWidth(2); + fitlpt->Draw("SAME"); + leg = new TLegend(0.5,0.38,0.85,0.53); + leg->SetBorderSize(0); + leg->SetFillColor(0); + leg->SetTextSize(0.05); + leg->SetTextFont(22); + leg->SetHeader("Lpt trigger pt cut"); + leg->AddEntry(fitlpt,"Ref.","l"); + leg->AddEntry(ptlpt,"New","l"); + leg->Draw("SAME"); + + c1->cd(3); + pthpt->SetMinimum(0.); + pthpt->SetMaximum(1.05); + pthpt->SetTitle(""); + pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); + pthpt->GetYaxis()->SetTitle("Efficiency"); + //pthpt->GetXaxis()->SetRange(3); + pthpt->Draw(""); + fithpt->SetParameters(0.627,0.393,0.855,0.0081); + fithpt->SetLineColor(2); + fithpt->SetLineWidth(2); + fithpt->Draw("SAME"); + leg = new TLegend(0.5,0.38,0.85,0.53); + leg->SetBorderSize(0); + leg->SetFillColor(0); + leg->SetTextSize(0.05); + leg->SetTextFont(22); + leg->SetHeader("Hpt trigger pt cut"); + leg->AddEntry(fithpt,"Ref.","l"); + leg->AddEntry(pthpt,"New","l"); + leg->Draw("SAME"); + + c1->SaveAs("MUONTriggerEfficiencyPt.gif"); + c1->SaveAs("MUONTriggerEfficiencyPt.eps"); + +} -- 2.43.0