/* $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
-*/
-
+/// \ingroup macros
+/// \file MUONTriggerEfficiencyPt.C
+/// \brief Macro to produce trigger single muon efficiency versus pt plots for the
+/// 2 pt cuts.
+///
+/// See full description on the \ref README_trigger page.
+///
+/// \author Fabien Guerin (LPC)
// ROOT includes
#include "TBranch.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLegend.h"
-
+#include "Riostream.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"
+#include "AliCDBManager.h"
// MUON includes
-#include "AliMUON.h"
-#include "AliMUONData.h"
+#include "AliMUONDataInterface.h"
+#include "AliMUONMCDataInterface.h"
+#include "AliMUONVHitStore.h"
+#include "AliMUONVTriggerStore.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)
return h0*TMath::TanH(h1)+par[3];
}
-void MUONTriggerEfficiencyPt(char filename[10]="galice.root")
+void MUONTriggerEfficiencyPt(const char* filenameSim="galice_sim.root",
+ const char* filenameRec="galice.root",
+ Bool_t readFromRP = 0)
{
// define style
st1->SetPadBottomMargin(0.15);
st1->cd();
- gROOT->ForceStyle();
+// 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;
+ Int_t nevents;
Double_t coincmuon=0;
- Double_t aptmuon=0;
Double_t lptmuon=0;
Double_t hptmuon=0;
- Double_t ptmu;
+ Double_t ptmu=0;
-// output file
- char digitdat[100];
- sprintf(digitdat,"MUONTriggerEfficiencyPt.out");
- FILE *fdat=fopen(digitdat,"w");
+ AliCDBManager* cdbManager = AliCDBManager::Instance();
+ cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ cdbManager->SetRun(0);
-// Creating Run Loader and openning file containing Hits
- AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
+ AliMUONMCDataInterface diSim(filenameSim);
+ AliMUONDataInterface diRec(filenameRec);
- if (RunLoader ==0x0) {
- printf(">>> Error : Error Opening %s file \n",filename);
+ if (!diSim.IsValid())
+ {
+ cerr << "Cannot access " << filenameSim << endl;
return;
}
- Int_t nevents = RunLoader->GetNumberOfEvents();
-
-// loading hits
- AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader");
- MUONLoader->LoadHits("READ");
- AliMUONData data_hits(MUONLoader,"MUON","MUON");
+ if (!diRec.IsValid())
+ {
+ cerr << "Cannot access " << filenameRec << endl;
+ return;
+ }
- TClonesArray * globalTrigger;
- AliMUONGlobalTrigger * gloTrg;
+ FILE* fdat = fopen("MUONTriggerEfficiencyPt.out","w");
+ if (!fdat)
+ {
+ cerr << "Cannot create output file" << endl;
+ return;
+ }
-// loading trigger
- MUONLoader->LoadDigits("READ");
- AliMUONData muondata(MUONLoader,"MUON","MUON");
+ nevents = diSim.NumberOfEvents();
+
+ AliRunLoader* runLoader = AliRunLoader::Open(filenameSim);
-// Loading kine
- RunLoader->LoadKinematics("READ");
-
for (Int_t ievent=0; ievent<nevents; ievent++) { // Event loop
- RunLoader->GetEvent(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; iparticle<nparticles; iparticle++) {
+
+ runLoader->GetEvent(ievent);
+ runLoader->LoadKinematics();
+ AliStack* stack = runLoader->Stack();
+
+ Int_t nparticles = (Int_t) stack->GetNtrack();
+
+ for (Int_t iparticle=0; iparticle<nparticles; iparticle++) {
particle = stack->Particle(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<nglobals; iglobal++) { // Global Trigger
- gloTrg = static_cast<AliMUONGlobalTrigger*>(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();
+// global trigger
+ TString tree("D");
+ if ( readFromRP ) tree = "R";
+
+ AliMUONVTriggerStore* triggerStore = diRec.TriggerStore(ievent,tree.Data());
+ AliMUONGlobalTrigger* gloTrg = triggerStore->Global();
+
+ if (gloTrg->SingleLpt()>=1 ) {
+ lptmuon++;
+ ptlpt->Fill(ptmu);
+ }
+ if (gloTrg->SingleHpt()>=1 ) {
+ hptmuon++;
+ pthpt->Fill(ptmu);
+ }
// 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();
-
+ ntracks = (Int_t) diSim.NumberOfTracks(ievent);
+
for (itrack=0; itrack<ntracks; itrack++) { // track loop
- data_hits.GetTrack(itrack);
+ AliMUONVHitStore* hitStore = diSim.HitStore(ievent,itrack);
+ AliMUONHit* mHit;
+ TIter next(hitStore->CreateIterator());
+
+ while ( ( mHit = static_cast<AliMUONHit*>(next()) ) )
+ {
+ Int_t Nch = mHit->Chamber();
+ Int_t hittrack = mHit->Track();
+ Float_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]++;
+ }
+ }
+ }
+ } // end track loop
- Int_t ihit, nhits;
- nhits = (Int_t) data_hits.Hits()->GetEntriesFast();
- AliMUONHit* mHit;
-
- for (ihit=0; ihit<nhits; ihit++) {
- mHit = static_cast<AliMUONHit*>(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");
-
+ } // end loop on event
+
+ if (coincmuon==0) {
+ cout << " >>> <E> coincmuon = 0 after event loop " << "\n";
+ cout << " >>> this probably means that input does not contain one (and only one) " << "\n";
+ cout << " >>> muon track per event as it should " << "\n";
+ cout << " >>> see README for further information " << "\n";
+ cout << " >>> exiting now ! " << "\n";
+ return;
+ }
+
+
+ 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," 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=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=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 *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(2,1);
+
+ c1->cd(1);
+ 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");
+ TLegend * leg = new TLegend(0.5,0.38,0.85,0.53);
+ 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(2);
+ 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");
+
}