Place the config and root file at the right place
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEefficiency.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
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 **************************************************************************/
15 //
16 // Task for Efficiency studies
17 // Used for testing classes AliHFEcontainer and AliHFEfilter
18 // Creates Efficiency Histograms
19 //
20 // Author:
21 //   Markus Fasel <M.Fasel@gsi.de>
22 //
23 #include <TAxis.h>
24 #include <TCanvas.h>
25 #include <TFile.h>
26 #include <TGraphErrors.h>
27 #include <TH1F.h>
28 #include <TH2.h>
29 #include <TIterator.h>
30 #include <TLegend.h>
31 #include <TObjArray.h>
32 #include <TPad.h>
33
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"
49
50 ClassImp(AliHFEefficiency)
51
52 AliHFEefficiency::AliHFEefficiency():
53   AliAnalysisTaskSE(),
54   fFilter(NULL),
55   fMCcut(NULL),
56   fAcceptanceCuts(NULL),
57   fEfficiency(NULL),
58   fOutput(NULL),
59   fCutTRD(kFALSE)
60 {
61   //
62   // Default constructor
63   //
64 }
65
66 AliHFEefficiency::AliHFEefficiency(const Char_t *name):
67   AliAnalysisTaskSE(name),
68   fFilter(NULL),
69   fMCcut(NULL),
70   fAcceptanceCuts(NULL),
71   fEfficiency(NULL),
72   fOutput(NULL),
73   fCutTRD(kFALSE)
74 {
75   //
76   // Main Constructor
77   //
78   DefineOutput(1, AliHFEcontainer::Class());
79   DefineOutput(2, TList::Class());
80
81   SetRunTerminate();
82 }
83
84 AliHFEefficiency::~AliHFEefficiency(){
85   //
86   // Destructor
87   //
88   if(fFilter) delete fFilter;
89   if(fEfficiency) delete fEfficiency;
90   if(fAcceptanceCuts) delete fAcceptanceCuts;
91   if(fOutput) delete fOutput;
92 }
93
94 void AliHFEefficiency::UserCreateOutputObjects(){
95   //
96   // Create the output objects
97   //
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);
109   // New container
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");
114
115   fFilter = new AliHFEtrackFilter("testfilter");
116   fFilter->GenerateCutSteps();
117   fFilter->InitCF(fEfficiency);
118   fMCcut = fFilter->GetMCSignalCuts();
119   if(fCutTRD){
120     AliHFEcutStep *hfeTRD = fFilter->GetCutStep("HFETRD");
121     AliHFEextraCuts *hfecuts = dynamic_cast<AliHFEextraCuts *>(hfeTRD->GetCut("HFETRDCuts"));
122     if(hfecuts) hfecuts->SetMinTrackletsTRD(4);
123   }
124   AliHFEcutStep *hfeITS = fFilter->GetCutStep("HFEITS");
125   AliHFEextraCuts *hfeitscuts = dynamic_cast<AliHFEextraCuts *>(hfeITS->GetCut("HFEPixelsCuts"));
126   if(hfeitscuts)hfeitscuts->SetRequireITSpixel(AliHFEextraCuts::kFirst);
127
128   AliHFEcutStep *primary = fFilter->GetCutStep("Primary");
129   AliCFTrackIsPrimaryCuts *primcuts = dynamic_cast<AliCFTrackIsPrimaryCuts *>(primary->GetCut("PrimaryCuts"));
130   if(primcuts){ 
131     primcuts->SetMaxDCAToVertexXY(3);
132     primcuts->SetMaxDCAToVertexZ(5);
133   }
134
135   fAcceptanceCuts = new AliCFAcceptanceCuts("Acceptance", "MC Acceptance Cuts");
136   fAcceptanceCuts->SetMinNHitITS(3);
137   fAcceptanceCuts->SetMinNHitTPC(2);
138   fAcceptanceCuts->SetMinNHitTRD(0);
139
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);
148 }
149
150 void AliHFEefficiency::UserExec(Option_t *){
151   //
152   // Event Loop
153   // Filter track, fill Efficiency container
154   //
155   AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fInputEvent);
156   if(!esdevent){
157     AliError("ESD Event required");
158     return;
159   }
160   fEfficiency->NewEvent();
161   fFilter->SetRecEvent(fInputEvent);
162   if(fMCEvent){
163     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
164    if ( ! mcH ) {
165      AliError("Cannot get MC truth event handler");
166      return;
167     }  
168     if(mcH &&(!mcH->InitOk())) return;
169     if(mcH &&(!mcH->TreeK())) return;
170     if(mcH &&(!mcH->TreeTR())) return;
171     fFilter->SetMC(fMCEvent);
172     FilterMC();
173   }
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);
184     else
185       fOutput->Fill("itspixel",0);
186     if(fMCEvent){
187       AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
188       if(mctrack){
189         Int_t motherLabel = mctrack->Particle()->GetFirstMother();
190         if(motherLabel){
191           AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherLabel));
192           if(mother)fOutput->Fill("mcmother", mother->Particle()->GetPdgCode());
193         }
194         fOutput->Fill("ptres", mctrack->Pt(), (track->Pt() - mctrack->Pt())/mctrack->Pt());
195       }
196     }
197   }
198   delete iter;
199   fFilter->Flush();
200   PostData(1, fEfficiency);
201   PostData(2, fOutput->GetList());
202 }
203
204 void AliHFEefficiency::Terminate(Option_t *){
205   //
206   // Evaluate Results
207   //
208   fEfficiency = dynamic_cast<AliHFEcontainer *>(GetOutputData(1));
209   if(!fEfficiency){
210     AliError("No Output data available");
211     return;
212   }
213
214   if(!IsRunningTerminate()) return;
215   PostProcess();
216
217   TList *l = dynamic_cast<TList *>(GetOutputData(2));
218   if(l) DrawPtResolution(l);
219 }
220
221 void AliHFEefficiency::FilterMC(){
222   //
223   // Cut MC tracks
224   //
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));
229     if(!track) continue;
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);
240   }
241 }
242
243 void AliHFEefficiency::Load(const char* filename){
244   //
245   // Load results for post processing
246   //
247   TFile *input = TFile::Open(filename);
248   AliHFEcontainer *ctmp = dynamic_cast<AliHFEcontainer *>(input->Get("Efficiency"));
249   if(ctmp) fEfficiency = dynamic_cast<AliHFEcontainer *>(ctmp->Clone());
250   input->Close();
251   delete input;
252 }
253
254 void AliHFEefficiency::PostProcess(){
255   //
256   // Calculate Efficiencies
257   //
258   fEfficiency->Print();
259
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);
265
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)));
269     ctmp->cd();
270     DrawSignalEfficiency(gridRec, crec, ivar);
271     
272     ctmp = new TCanvas(Form("cCutEff%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s)", crec->GetVarTitle(ivar)));
273     ctmp->cd();
274     DrawCutEfficiency(gridRec, crec, ivar);
275   
276     ctmp = new TCanvas(Form("cEffSignalMC%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
277     ctmp->cd();
278     DrawSignalEfficiency(gridMC, cmc, ivar);
279
280     ctmp = new TCanvas(Form("cCutEffMC%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
281     ctmp->cd();
282     DrawCutEfficiency(gridMC, cmc, ivar);
283   } 
284
285   CalculatePTsmearing();
286
287   delete gridRec;
288   delete gridMC;
289   delete crec;
290   delete cmc;
291 }
292
293 void AliHFEefficiency::DrawSignalEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
294   //
295   // Efficiency of a cut step with respect to the signal
296   //
297   TH1 *hEff = NULL;
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");
310     kFirst = kFALSE;
311     leg->AddEntry(hEff, cont->GetStepTitle(istep), "p");
312   }
313   leg->SetBorderSize(0);
314   leg->SetFillStyle(0);
315   leg->Draw();
316   gPad->Update();
317 }
318
319 void AliHFEefficiency::DrawCutEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
320   //
321   // Efficiency of a cut step with respect to the previous one
322   //
323   Bool_t kFirst = kTRUE;
324   TH1 *hEff = NULL;
325   TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
326   TString effTitle;
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");
342     kFirst = kFALSE;
343     leg->AddEntry(hEff, cont->GetStepTitle(istep), "p");
344   }
345   leg->SetBorderSize(0);
346   leg->SetFillStyle(0);
347   leg->Draw();
348   gPad->Update();
349 }
350
351 void AliHFEefficiency::CalculatePTsmearing(){
352   //
353   // Efficiency of a cut step with respect to the same cut step filled with MC 
354   // Information
355   //
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;
359   TH1 *hEff = NULL;
360   TCanvas *c2 = new TCanvas("cSmear", "pT modification");
361   c2->cd();
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");
371     kFirst = kFALSE;
372     leg1->AddEntry(hEff, cont->GetStepTitle(istep), "p");
373   }
374   leg1->Draw();
375   gPad->Update();
376
377   delete cont;
378   delete grid;
379 }
380
381 void AliHFEefficiency::DrawPtResolution(const TList * const l){
382   //
383   // Draw pt resolution
384   //
385   TH2 *h2 = dynamic_cast<TH2 *>(l->FindObject("ptres"));
386   if(!h2) return;
387   TGraphErrors *mean = new TGraphErrors(h2->GetNbinsX());
388   TGraphErrors *rms = new TGraphErrors(h2->GetNbinsX());
389
390   Double_t pt = 0., pterr = 0., mymean = 0., myrms = 0., mymeanerr = 0., myrmserr = 0.;
391   TH1 *htmp = NULL;
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();
400     delete htmp;
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);
405   }
406
407   TCanvas *c1 = new TCanvas("cPtShift", "pT Shift");
408   c1->cd();
409   mean->SetMarkerStyle(22);
410   mean->SetMarkerColor(kBlue);
411   mean->Draw("ape");
412   rms->SetMarkerStyle(22);
413   rms->SetMarkerColor(kBlack);
414   rms->Draw("pesame");
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");
418   leg->Draw();
419   gPad->Update();
420
421