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());
164 if(mcH &&(!mcH->InitOk())) return;
165 if(mcH &&(!mcH->TreeK())) return;
166 if(mcH &&(!mcH->TreeTR())) return;
167 fFilter->SetMC(fMCEvent);
170 fFilter->FilterTracks(esdevent);
171 TObjArray *tracks = fFilter->GetFilteredTracks();
172 TIterator *iter = tracks->MakeIterator();
173 fOutput->Fill("hNtracks", tracks->GetEntriesFast());
174 AliESDtrack *track = NULL;
175 while((track = dynamic_cast<AliESDtrack *>(iter->Next()))){
176 fOutput->Fill("hPt", track->Pt());
177 fOutput->Fill("kinkIndex", track->GetKinkIndex(0));
178 if(TESTBIT(track->GetITSClusterMap(), 0))
179 fOutput->Fill("itspixel",1);
181 fOutput->Fill("itspixel",0);
183 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
185 Int_t motherLabel = mctrack->Particle()->GetFirstMother();
187 AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherLabel));
188 if(mother)fOutput->Fill("mcmother", mother->Particle()->GetPdgCode());
190 fOutput->Fill("ptres", mctrack->Pt(), (track->Pt() - mctrack->Pt())/mctrack->Pt());
196 PostData(1, fEfficiency);
197 PostData(2, fOutput->GetList());
200 void AliHFEefficiency::Terminate(Option_t *){
204 fEfficiency = dynamic_cast<AliHFEcontainer *>(GetOutputData(1));
206 AliError("No Output data available");
210 if(!IsRunningTerminate()) return;
213 TList *l = dynamic_cast<TList *>(GetOutputData(2));
214 if(l) DrawPtResolution(l);
217 void AliHFEefficiency::FilterMC(){
221 Double_t cont[4] = {0., 0., 0., 0.};
222 AliMCParticle *track = NULL;
223 for(Int_t itrack = 0; itrack < fMCEvent->GetNumberOfTracks(); itrack++){
224 track = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrack));
226 if(!fMCcut->IsSelected(track)) continue;
227 cont[0] = track->Pt();
228 cont[1] = track->Eta();
229 cont[2] = track->Phi();
230 cont[3] = track->Charge()/3;
231 //AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(track->Particle()->GetFirstMother()));
232 //if(TMath::Abs(mother->Particle()->GetPdgCode()) != 443) continue;
233 fEfficiency->FillCFContainer("MCFilter", 0, cont);
234 if(!fAcceptanceCuts->IsSelected(track)) continue;
235 fEfficiency->FillCFContainer("MCFilter", 1, cont);
239 void AliHFEefficiency::Load(const char* filename){
241 // Load results for post processing
243 TFile *input = TFile::Open(filename);
244 AliHFEcontainer *cin = dynamic_cast<AliHFEcontainer *>(input->Get("Efficiency"));
245 if(cin) fEfficiency = dynamic_cast<AliHFEcontainer *>(cin->Clone());
250 void AliHFEefficiency::PostProcess(){
252 // Calculate Efficiencies
254 fEfficiency->Print();
256 // Prepare containers
257 AliCFContainer *crec = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signal");
258 AliCFContainer *cmc = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signalMC");
259 AliCFEffGrid *gridRec = new AliCFEffGrid("effrec", "Efficiency Calculation using rec vars", *crec);
260 AliCFEffGrid *gridMC = new AliCFEffGrid("effmc", "Efficiency Calculation using mc vars", *cmc);
262 TCanvas *ctmp = NULL;
263 for(Int_t ivar = 0; ivar < crec->GetNVar(); ivar++){
264 ctmp = new TCanvas(Form("cEffSignal%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s)", crec->GetVarTitle(ivar)));
266 DrawSignalEfficiency(gridRec, crec, ivar);
268 ctmp = new TCanvas(Form("cCutEff%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s)", crec->GetVarTitle(ivar)));
270 DrawCutEfficiency(gridRec, crec, ivar);
272 ctmp = new TCanvas(Form("cEffSignalMC%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
274 DrawSignalEfficiency(gridMC, cmc, ivar);
276 ctmp = new TCanvas(Form("cCutEffMC%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
278 DrawCutEfficiency(gridMC, cmc, ivar);
281 CalculatePTsmearing();
289 void AliHFEefficiency::DrawSignalEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
291 // Efficiency of a cut step with respect to the signal
294 Bool_t kFirst = kTRUE;
295 TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
296 for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
297 eff->CalculateEfficiency(istep, 0);
298 hEff = eff->Project(var);
299 hEff->SetTitle(Form("Signal Efficiency (%s)", cont->GetVarTitle(var)));
300 hEff->GetYaxis()->SetRangeUser(0., 1.6);
301 hEff->GetYaxis()->SetTitle("Efficiency");
302 hEff->SetStats(kFALSE);
303 hEff->SetMarkerStyle(20+istep);
304 hEff->SetMarkerColor(kBlue - 5 + istep);
305 hEff->Draw(kFirst ? "ep" : "epsame");
307 leg->AddEntry(hEff, cont->GetStepTitle(istep), "p");
309 leg->SetBorderSize(0);
310 leg->SetFillStyle(0);
315 void AliHFEefficiency::DrawCutEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
317 // Efficiency of a cut step with respect to the previous one
319 Bool_t kFirst = kTRUE;
321 TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
323 Int_t start = 0, length = 0;
324 for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
325 eff->CalculateEfficiency(istep, istep -1);
326 hEff = eff->Project(var);
327 effTitle = hEff->GetTitle();
328 start = effTitle.Index("projected"); length = effTitle.Length() - start;
329 effTitle.Remove(start - 1, length + 1);
330 effTitle.ReplaceAll("Efficiency: ", "");
331 hEff->SetTitle(Form("Cut Efficiency (%s)", cont->GetVarTitle(var)));
332 hEff->GetYaxis()->SetRangeUser(0., 1.6);
333 hEff->GetYaxis()->SetTitle("Efficiency");
334 hEff->SetStats(kFALSE);
335 hEff->SetMarkerStyle(20+istep);
336 hEff->SetMarkerColor(kBlue - 5 + istep);
337 hEff->Draw(kFirst ? "ep" : "epsame");
339 leg->AddEntry(hEff, effTitle.Data(), "p");
341 leg->SetBorderSize(0);
342 leg->SetFillStyle(0);
347 void AliHFEefficiency::CalculatePTsmearing(){
349 // Efficiency of a cut step with respect to the same cut step filled with MC
352 AliCFContainer *cont = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "testfilter_container_signal:testfilter_container_signalMC");
353 AliCFEffGrid *grid = new AliCFEffGrid("effc", "Efficiency Calculation", *cont);
354 Bool_t kFirst = kTRUE;
356 TCanvas *c2 = new TCanvas("cSmear", "pT modification");
358 TLegend *leg1 = new TLegend(0.5, 0.7, 0.89, 0.89);
359 Int_t nStep = cont->GetNStep()/2;
360 for(Int_t istep = 0; istep < nStep; istep++){
361 grid->CalculateEfficiency(istep, istep + nStep);
362 hEff = grid->Project(0);
363 hEff->SetStats(kFALSE);
364 hEff->SetMarkerStyle(20+istep);
365 hEff->SetMarkerColor(kBlue - 5 + istep);
366 hEff->Draw(kFirst ? "ep" : "epsame");
368 leg1->AddEntry(hEff, cont->GetStepTitle(istep), "p");
377 void AliHFEefficiency::DrawPtResolution(const TList * const l){
379 // Draw pt resolution
381 TH2 *h2 = dynamic_cast<TH2 *>(l->FindObject("ptres"));
383 TGraphErrors *mean = new TGraphErrors(h2->GetNbinsX());
384 TGraphErrors *rms = new TGraphErrors(h2->GetNbinsX());
386 Double_t pt = 0., pterr = 0., mymean = 0., myrms = 0., mymeanerr = 0., myrmserr = 0.;
388 for(Int_t imom = 1; imom <= h2->GetXaxis()->GetLast(); imom++){
389 pt = h2->GetXaxis()->GetBinCenter(imom);
390 pterr = h2->GetXaxis()->GetBinWidth(imom)/2.;
391 htmp = h2->ProjectionY("py", imom, imom);
392 mymean = htmp->GetMean();
393 myrms = htmp->GetRMS();
394 mymeanerr = htmp->GetMeanError();
395 myrmserr = htmp->GetRMSError();
397 mean->SetPoint(imom -1, pt, mymean);
398 rms->SetPoint(imom -1, pt, myrms);
399 mean->SetPointError(imom-1,pterr,mymeanerr);
400 rms->SetPointError(imom-1,pterr,myrmserr);
403 TCanvas *c1 = new TCanvas("cPtShift", "pT Shift");
405 mean->SetMarkerStyle(22);
406 mean->SetMarkerColor(kBlue);
408 rms->SetMarkerStyle(22);
409 rms->SetMarkerColor(kBlack);
411 TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
412 leg->AddEntry(mean, "Mean", "lp");
413 leg->AddEntry(rms, "RMS", "lp");