1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Task for Efficiency studies
17 // Used for testing classes AliHFEcontainer and AliHFEfilter
18 // Creates Efficiency Histograms
21 // Markus Fasel <M.Fasel@gsi.de>
26 #include <TGraphErrors.h>
29 #include <TIterator.h>
31 #include <TObjArray.h>
34 #include "AliAnalysisManager.h"
35 #include "AliCFAcceptanceCuts.h"
36 #include "AliCFTrackIsPrimaryCuts.h"
37 #include "AliCFContainer.h"
38 #include "AliCFEffGrid.h"
39 #include "AliESDEvent.h"
40 #include "AliHFEcollection.h"
41 #include "AliHFEcontainer.h"
42 #include "AliHFEcutStep.h"
43 #include "AliHFEefficiency.h"
44 #include "AliHFEextraCuts.h"
45 #include "AliHFEtrackFilter.h"
46 #include "AliHFEtools.h"
47 #include "AliMCEvent.h"
48 #include "AliMCEventHandler.h"
50 ClassImp(AliHFEefficiency)
52 AliHFEefficiency::AliHFEefficiency():
56 fAcceptanceCuts(NULL),
62 // Default constructor
66 AliHFEefficiency::AliHFEefficiency(const Char_t *name):
67 AliAnalysisTaskSE(name),
70 fAcceptanceCuts(NULL),
78 DefineOutput(1, AliHFEcontainer::Class());
79 DefineOutput(2, TList::Class());
84 AliHFEefficiency::~AliHFEefficiency(){
88 if(fFilter) delete fFilter;
89 if(fEfficiency) delete fEfficiency;
90 if(fAcceptanceCuts) delete fAcceptanceCuts;
91 if(fOutput) delete fOutput;
94 void AliHFEefficiency::UserCreateOutputObjects(){
96 // Create the output objects
98 fEfficiency = new AliHFEcontainer("Efficiency");
99 fEfficiency->SetNumberOfVariables(4);
100 // InitializeVariables
101 fEfficiency->SetVariableName(0, "pt");
102 fEfficiency->MakeLogarithmicBinning(0, 40, 0.1, 10);
103 fEfficiency->SetVariableName(1, "Eta");
104 fEfficiency->MakeLinearBinning(1, 20, -0.9, 0.9);
105 fEfficiency->SetVariableName(2, "phi");
106 fEfficiency->MakeLinearBinning(2 , 18, 0, 2 * TMath::Pi());
107 fEfficiency->SetVariableName(3, "Charge");
108 fEfficiency->MakeLinearBinning(3, 2, -1.1, 1.1);
110 fEfficiency->CreateContainer("MCFilter", "Efficiency for MC Studies", 2);
111 AliCFContainer *cont = fEfficiency->GetCFContainer("MCFilter");
112 cont->SetStepTitle(0, "MC Signal");
113 cont->SetStepTitle(1, "MC Signal in acceptance");
115 fFilter = new AliHFEtrackFilter("testfilter");
116 fFilter->GenerateCutSteps();
117 fFilter->InitCF(fEfficiency);
118 fMCcut = fFilter->GetMCSignalCuts();
120 AliHFEcutStep *hfeTRD = fFilter->GetCutStep("HFETRD");
121 AliHFEextraCuts *hfecuts = dynamic_cast<AliHFEextraCuts *>(hfeTRD->GetCut("HFETRDCuts"));
122 if(hfecuts) hfecuts->SetMinTrackletsTRD(4);
124 AliHFEcutStep *hfeITS = fFilter->GetCutStep("HFEITS");
125 AliHFEextraCuts *hfeitscuts = dynamic_cast<AliHFEextraCuts *>(hfeITS->GetCut("HFEPixelsCuts"));
126 if(hfeitscuts)hfeitscuts->SetRequireITSpixel(AliHFEextraCuts::kFirst);
128 AliHFEcutStep *primary = fFilter->GetCutStep("Primary");
129 AliCFTrackIsPrimaryCuts *primcuts = dynamic_cast<AliCFTrackIsPrimaryCuts *>(primary->GetCut("PrimaryCuts"));
131 primcuts->SetMaxDCAToVertexXY(3);
132 primcuts->SetMaxDCAToVertexZ(5);
135 fAcceptanceCuts = new AliCFAcceptanceCuts("Acceptance", "MC Acceptance Cuts");
136 fAcceptanceCuts->SetMinNHitITS(3);
137 fAcceptanceCuts->SetMinNHitTPC(2);
138 fAcceptanceCuts->SetMinNHitTRD(0);
140 // create additional histos for testing purpose
141 fOutput = new AliHFEcollection("histos", "QA histograms");
142 fOutput->CreateTH1F("hNtracks","Number of Tracks per Event", 100, 0, 100);
143 fOutput->CreateTH1F("hPt","Pt of the tracks", 40, 0.1, 10, 0);
144 fOutput->CreateTH1F("kinkIndex", "Kink Index", 400, -200, 200);
145 fOutput->CreateTH1F("itspixel", "ITS PIXEL", 2, 0, 2);
146 fOutput->CreateTH1F("mcmother", "Mother PDG", 1000, -500, 500);
147 fOutput->CreateTH2F("ptres", "Pt Resolution", 40, 0.1, 10, 200, -1.5, 1.5, 0);
150 void AliHFEefficiency::UserExec(Option_t *){
153 // Filter track, fill Efficiency container
155 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fInputEvent);
157 AliError("ESD Event required");
160 fEfficiency->NewEvent();
161 fFilter->SetRecEvent(fInputEvent);
163 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
165 AliError("Cannot get MC truth event handler");
168 if(mcH &&(!mcH->InitOk())) return;
169 if(mcH &&(!mcH->TreeK())) return;
170 if(mcH &&(!mcH->TreeTR())) return;
171 fFilter->SetMC(fMCEvent);
174 fFilter->FilterTracks(esdevent);
175 TObjArray *tracks = fFilter->GetFilteredTracks();
176 TIterator *iter = tracks->MakeIterator();
177 fOutput->Fill("hNtracks", tracks->GetEntriesFast());
178 AliESDtrack *track = NULL;
179 while((track = dynamic_cast<AliESDtrack *>(iter->Next()))){
180 fOutput->Fill("hPt", track->Pt());
181 fOutput->Fill("kinkIndex", track->GetKinkIndex(0));
182 if(TESTBIT(track->GetITSClusterMap(), 0))
183 fOutput->Fill("itspixel",1);
185 fOutput->Fill("itspixel",0);
187 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
189 Int_t motherLabel = mctrack->Particle()->GetFirstMother();
191 AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherLabel));
192 if(mother)fOutput->Fill("mcmother", mother->Particle()->GetPdgCode());
194 fOutput->Fill("ptres", mctrack->Pt(), (track->Pt() - mctrack->Pt())/mctrack->Pt());
200 PostData(1, fEfficiency);
201 PostData(2, fOutput->GetList());
204 void AliHFEefficiency::Terminate(Option_t *){
208 fEfficiency = dynamic_cast<AliHFEcontainer *>(GetOutputData(1));
210 AliError("No Output data available");
214 if(!IsRunningTerminate()) return;
217 TList *l = dynamic_cast<TList *>(GetOutputData(2));
218 if(l) DrawPtResolution(l);
221 void AliHFEefficiency::FilterMC(){
225 Double_t cont[4] = {0., 0., 0., 0.};
226 AliMCParticle *track = NULL;
227 for(Int_t itrack = 0; itrack < fMCEvent->GetNumberOfTracks(); itrack++){
228 track = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrack));
230 if(!fMCcut->IsSelected(track)) continue;
231 cont[0] = track->Pt();
232 cont[1] = track->Eta();
233 cont[2] = track->Phi();
234 cont[3] = track->Charge()/3;
235 //AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(track->Particle()->GetFirstMother()));
236 //if(TMath::Abs(mother->Particle()->GetPdgCode()) != 443) continue;
237 fEfficiency->FillCFContainer("MCFilter", 0, cont);
238 if(!fAcceptanceCuts->IsSelected(track)) continue;
239 fEfficiency->FillCFContainer("MCFilter", 1, cont);
243 void AliHFEefficiency::Load(const char* filename){
245 // Load results for post processing
247 TFile *input = TFile::Open(filename);
248 AliHFEcontainer *ctmp = dynamic_cast<AliHFEcontainer *>(input->Get("Efficiency"));
249 if(ctmp) fEfficiency = dynamic_cast<AliHFEcontainer *>(ctmp->Clone());
254 void AliHFEefficiency::PostProcess(){
256 // Calculate Efficiencies
258 fEfficiency->Print();
260 // Prepare containers
261 AliCFContainer *crec = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signal");
262 AliCFContainer *cmc = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signalMC");
263 AliCFEffGrid *gridRec = new AliCFEffGrid("effrec", "Efficiency Calculation using rec vars", *crec);
264 AliCFEffGrid *gridMC = new AliCFEffGrid("effmc", "Efficiency Calculation using mc vars", *cmc);
266 TCanvas *ctmp = NULL;
267 for(Int_t ivar = 0; ivar < crec->GetNVar(); ivar++){
268 ctmp = new TCanvas(Form("cEffSignal%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s)", crec->GetVarTitle(ivar)));
270 DrawSignalEfficiency(gridRec, crec, ivar);
272 ctmp = new TCanvas(Form("cCutEff%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s)", crec->GetVarTitle(ivar)));
274 DrawCutEfficiency(gridRec, crec, ivar);
276 ctmp = new TCanvas(Form("cEffSignalMC%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
278 DrawSignalEfficiency(gridMC, cmc, ivar);
280 ctmp = new TCanvas(Form("cCutEffMC%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
282 DrawCutEfficiency(gridMC, cmc, ivar);
285 CalculatePTsmearing();
293 void AliHFEefficiency::DrawSignalEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
295 // Efficiency of a cut step with respect to the signal
298 Bool_t kFirst = kTRUE;
299 TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
300 for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
301 eff->CalculateEfficiency(istep, 0);
302 hEff = eff->Project(var);
303 hEff->SetTitle(Form("Signal Efficiency (%s)", cont->GetVarTitle(var)));
304 hEff->GetYaxis()->SetRangeUser(0., 1.6);
305 hEff->GetYaxis()->SetTitle("Efficiency");
306 hEff->SetStats(kFALSE);
307 hEff->SetMarkerStyle(20+istep);
308 hEff->SetMarkerColor(kBlue - 5 + istep);
309 hEff->Draw(kFirst ? "ep" : "epsame");
311 leg->AddEntry(hEff, cont->GetStepTitle(istep), "p");
313 leg->SetBorderSize(0);
314 leg->SetFillStyle(0);
319 void AliHFEefficiency::DrawCutEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
321 // Efficiency of a cut step with respect to the previous one
323 Bool_t kFirst = kTRUE;
325 TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
327 Int_t start = 0, length = 0;
328 for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
329 eff->CalculateEfficiency(istep, istep -1);
330 hEff = eff->Project(var);
331 effTitle = hEff->GetTitle();
332 start = effTitle.Index("projected"); length = effTitle.Length() - start;
333 effTitle.Remove(start - 1, length + 1);
334 effTitle.ReplaceAll("Efficiency: ", "");
335 hEff->SetTitle(Form("Cut Efficiency (%s)", cont->GetVarTitle(var)));
336 hEff->GetYaxis()->SetRangeUser(0., 1.6);
337 hEff->GetYaxis()->SetTitle("Efficiency");
338 hEff->SetStats(kFALSE);
339 hEff->SetMarkerStyle(20+istep);
340 hEff->SetMarkerColor(kBlue - 5 + istep);
341 hEff->Draw(kFirst ? "ep" : "epsame");
343 leg->AddEntry(hEff, cont->GetStepTitle(istep), "p");
345 leg->SetBorderSize(0);
346 leg->SetFillStyle(0);
351 void AliHFEefficiency::CalculatePTsmearing(){
353 // Efficiency of a cut step with respect to the same cut step filled with MC
356 AliCFContainer *cont = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "testfilter_container_signal:testfilter_container_signalMC");
357 AliCFEffGrid *grid = new AliCFEffGrid("effc", "Efficiency Calculation", *cont);
358 Bool_t kFirst = kTRUE;
360 TCanvas *c2 = new TCanvas("cSmear", "pT modification");
362 TLegend *leg1 = new TLegend(0.5, 0.7, 0.89, 0.89);
363 Int_t nStep = cont->GetNStep()/2;
364 for(Int_t istep = 0; istep < nStep; istep++){
365 grid->CalculateEfficiency(istep, istep + nStep);
366 hEff = grid->Project(0);
367 hEff->SetStats(kFALSE);
368 hEff->SetMarkerStyle(20+istep);
369 hEff->SetMarkerColor(kBlue - 5 + istep);
370 hEff->Draw(kFirst ? "ep" : "epsame");
372 leg1->AddEntry(hEff, cont->GetStepTitle(istep), "p");
381 void AliHFEefficiency::DrawPtResolution(const TList * const l){
383 // Draw pt resolution
385 TH2 *h2 = dynamic_cast<TH2 *>(l->FindObject("ptres"));
387 TGraphErrors *mean = new TGraphErrors(h2->GetNbinsX());
388 TGraphErrors *rms = new TGraphErrors(h2->GetNbinsX());
390 Double_t pt = 0., pterr = 0., mymean = 0., myrms = 0., mymeanerr = 0., myrmserr = 0.;
392 for(Int_t imom = 1; imom <= h2->GetXaxis()->GetLast(); imom++){
393 pt = h2->GetXaxis()->GetBinCenter(imom);
394 pterr = h2->GetXaxis()->GetBinWidth(imom)/2.;
395 htmp = h2->ProjectionY("py", imom, imom);
396 mymean = htmp->GetMean();
397 myrms = htmp->GetRMS();
398 mymeanerr = htmp->GetMeanError();
399 myrmserr = htmp->GetRMSError();
401 mean->SetPoint(imom -1, pt, mymean);
402 rms->SetPoint(imom -1, pt, myrms);
403 mean->SetPointError(imom-1,pterr,mymeanerr);
404 rms->SetPointError(imom-1,pterr,myrmserr);
407 TCanvas *c1 = new TCanvas("cPtShift", "pT Shift");
409 mean->SetMarkerStyle(22);
410 mean->SetMarkerColor(kBlue);
412 rms->SetMarkerStyle(22);
413 rms->SetMarkerColor(kBlack);
415 TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
416 leg->AddEntry(mean, "Mean", "lp");
417 leg->AddEntry(rms, "RMS", "lp");