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 fEfficiency->NewEvent();
156 fFilter->SetRecEvent(fInputEvent);
158 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
159 if(!mcH->InitOk()) return;
160 if(!mcH->TreeK()) return;
161 if(!mcH->TreeTR()) return;
162 fFilter->SetMC(fMCEvent);
165 fFilter->FilterTracks(dynamic_cast<AliESDEvent *>(fInputEvent));
166 TObjArray *tracks = fFilter->GetFilteredTracks();
167 TIterator *iter = tracks->MakeIterator();
168 fOutput->Fill("hNtracks", tracks->GetEntriesFast());
169 AliESDtrack *track = NULL;
170 while((track = dynamic_cast<AliESDtrack *>(iter->Next()))){
171 fOutput->Fill("hPt", track->Pt());
172 fOutput->Fill("kinkIndex", track->GetKinkIndex(0));
173 if(TESTBIT(track->GetITSClusterMap(), 0))
174 fOutput->Fill("itspixel",1);
176 fOutput->Fill("itspixel",0);
178 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
180 Int_t motherLabel = mctrack->Particle()->GetFirstMother();
182 AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherLabel));
183 if(mother)fOutput->Fill("mcmother", mother->Particle()->GetPdgCode());
185 fOutput->Fill("ptres", mctrack->Pt(), (track->Pt() - mctrack->Pt())/mctrack->Pt());
191 PostData(1, fEfficiency);
192 PostData(2, fOutput->GetList());
195 void AliHFEefficiency::Terminate(Option_t *){
199 fEfficiency = dynamic_cast<AliHFEcontainer *>(GetOutputData(1));
201 AliError("No Output data available");
205 if(!IsRunningTerminate()) return;
208 TList *l = dynamic_cast<TList *>(GetOutputData(2));
209 if(l) DrawPtResolution(l);
212 void AliHFEefficiency::FilterMC(){
216 Double_t cont[4] = {0., 0., 0., 0.};
217 AliMCParticle *track = NULL;
218 for(Int_t itrack = 0; itrack < fMCEvent->GetNumberOfTracks(); itrack++){
219 track = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrack));
220 if(!fMCcut->IsSelected(track)) continue;
221 cont[0] = track->Pt();
222 cont[1] = track->Eta();
223 cont[2] = track->Phi();
224 cont[3] = track->Charge()/3;
225 //AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(track->Particle()->GetFirstMother()));
226 //if(TMath::Abs(mother->Particle()->GetPdgCode()) != 443) continue;
227 fEfficiency->FillCFContainer("MCFilter", 0, cont);
228 if(!fAcceptanceCuts->IsSelected(track)) continue;
229 fEfficiency->FillCFContainer("MCFilter", 1, cont);
233 void AliHFEefficiency::Load(const char* filename){
235 // Load results for post processing
237 TFile *input = TFile::Open(filename);
238 AliHFEcontainer *cin = dynamic_cast<AliHFEcontainer *>(input->Get("Efficiency"));
239 if(cin) fEfficiency = dynamic_cast<AliHFEcontainer *>(cin->Clone());
244 void AliHFEefficiency::PostProcess(){
246 // Calculate Efficiencies
248 fEfficiency->Print();
250 // Prepare containers
251 AliCFContainer *crec = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signal");
252 AliCFContainer *cmc = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signalMC");
253 AliCFEffGrid *gridRec = new AliCFEffGrid("effrec", "Efficiency Calculation using rec vars", *crec);
254 AliCFEffGrid *gridMC = new AliCFEffGrid("effmc", "Efficiency Calculation using mc vars", *cmc);
256 TCanvas *ctmp = NULL;
257 for(Int_t ivar = 0; ivar < crec->GetNVar(); ivar++){
258 ctmp = new TCanvas(Form("cEffSignal%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s)", crec->GetVarTitle(ivar)));
260 DrawSignalEfficiency(gridRec, crec, ivar);
262 ctmp = new TCanvas(Form("cCutEff%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s)", crec->GetVarTitle(ivar)));
264 DrawCutEfficiency(gridRec, crec, ivar);
266 ctmp = new TCanvas(Form("cEffSignalMC%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
268 DrawSignalEfficiency(gridMC, cmc, ivar);
270 ctmp = new TCanvas(Form("cCutEffMC%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
272 DrawCutEfficiency(gridMC, cmc, ivar);
275 CalculatePTsmearing();
283 void AliHFEefficiency::DrawSignalEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
285 // Efficiency of a cut step with respect to the signal
288 Bool_t kFirst = kTRUE;
289 TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
290 for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
291 eff->CalculateEfficiency(istep, 0);
292 hEff = eff->Project(var);
293 hEff->SetTitle(Form("Signal Efficiency (%s)", cont->GetVarTitle(var)));
294 hEff->GetYaxis()->SetRangeUser(0., 1.6);
295 hEff->GetYaxis()->SetTitle("Efficiency");
296 hEff->SetStats(kFALSE);
297 hEff->SetMarkerStyle(20+istep);
298 hEff->SetMarkerColor(kBlue - 5 + istep);
299 hEff->Draw(kFirst ? "ep" : "epsame");
301 leg->AddEntry(hEff, cont->GetStepTitle(istep), "p");
303 leg->SetBorderSize(0);
304 leg->SetFillStyle(0);
309 void AliHFEefficiency::DrawCutEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
311 // Efficiency of a cut step with respect to the previous one
313 Bool_t kFirst = kTRUE;
315 TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
317 Int_t start = 0, length = 0;
318 for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
319 eff->CalculateEfficiency(istep, istep -1);
320 hEff = eff->Project(var);
321 effTitle = hEff->GetTitle();
322 start = effTitle.Index("projected"); length = effTitle.Length() - start;
323 effTitle.Remove(start - 1, length + 1);
324 effTitle.ReplaceAll("Efficiency: ", "");
325 hEff->SetTitle(Form("Cut Efficiency (%s)", cont->GetVarTitle(var)));
326 hEff->GetYaxis()->SetRangeUser(0., 1.6);
327 hEff->GetYaxis()->SetTitle("Efficiency");
328 hEff->SetStats(kFALSE);
329 hEff->SetMarkerStyle(20+istep);
330 hEff->SetMarkerColor(kBlue - 5 + istep);
331 hEff->Draw(kFirst ? "ep" : "epsame");
333 leg->AddEntry(hEff, effTitle.Data(), "p");
335 leg->SetBorderSize(0);
336 leg->SetFillStyle(0);
341 void AliHFEefficiency::CalculatePTsmearing(){
343 // Efficiency of a cut step with respect to the same cut step filled with MC
346 AliCFContainer *cont = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "testfilter_container_signal:testfilter_container_signalMC");
347 AliCFEffGrid *grid = new AliCFEffGrid("effc", "Efficiency Calculation", *cont);
348 Bool_t kFirst = kTRUE;
350 TCanvas *c2 = new TCanvas("cSmear", "pT modification");
352 TLegend *leg1 = new TLegend(0.5, 0.7, 0.89, 0.89);
353 Int_t nStep = cont->GetNStep()/2;
354 for(Int_t istep = 0; istep < nStep; istep++){
355 grid->CalculateEfficiency(istep, istep + nStep);
356 hEff = grid->Project(0);
357 hEff->SetStats(kFALSE);
358 hEff->SetMarkerStyle(20+istep);
359 hEff->SetMarkerColor(kBlue - 5 + istep);
360 hEff->Draw(kFirst ? "ep" : "epsame");
362 leg1->AddEntry(hEff, cont->GetStepTitle(istep), "p");
371 void AliHFEefficiency::DrawPtResolution(const TList * const l){
373 // Draw pt resolution
375 TH2 *h2 = dynamic_cast<TH2 *>(l->FindObject("ptres"));
377 TGraphErrors *mean = new TGraphErrors(h2->GetNbinsX());
378 TGraphErrors *rms = new TGraphErrors(h2->GetNbinsX());
380 Double_t pt = 0., pterr = 0., mymean = 0., myrms = 0., mymeanerr = 0., myrmserr = 0.;
382 for(Int_t imom = 1; imom <= h2->GetXaxis()->GetLast(); imom++){
383 pt = h2->GetXaxis()->GetBinCenter(imom);
384 pterr = h2->GetXaxis()->GetBinWidth(imom)/2.;
385 htmp = h2->ProjectionY("py", imom, imom);
386 mymean = htmp->GetMean();
387 myrms = htmp->GetRMS();
388 mymeanerr = htmp->GetMeanError();
389 myrmserr = htmp->GetRMSError();
391 mean->SetPoint(imom -1, pt, mymean);
392 rms->SetPoint(imom -1, pt, myrms);
393 mean->SetPointError(imom-1,pterr,mymeanerr);
394 rms->SetPointError(imom-1,pterr,myrmserr);
397 TCanvas *c1 = new TCanvas("cPtShift", "pT Shift");
399 mean->SetMarkerStyle(22);
400 mean->SetMarkerColor(kBlue);
402 rms->SetMarkerStyle(22);
403 rms->SetMarkerColor(kBlack);
405 TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
406 leg->AddEntry(mean, "Mean", "lp");
407 leg->AddEntry(rms, "RMS", "lp");