Major update of the HFE package (comments inside the code
[u/mrichter/AliRoot.git] / PWG3 / 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     hfecuts->SetMinTrackletsTRD(4);
123   }
124   AliHFEcutStep *hfeITS = fFilter->GetCutStep("HFEITS");
125   AliHFEextraCuts *hfeitscuts = dynamic_cast<AliHFEextraCuts *>(hfeITS->GetCut("HFEPixelsCuts"));
126   hfeitscuts->SetRequireITSpixel(AliHFEextraCuts::kFirst);
127
128   AliHFEcutStep *primary = fFilter->GetCutStep("Primary");
129   AliCFTrackIsPrimaryCuts *primcuts = dynamic_cast<AliCFTrackIsPrimaryCuts *>(primary->GetCut("PrimaryCuts"));
130   primcuts->SetMaxDCAToVertexXY(3);
131   primcuts->SetMaxDCAToVertexZ(5);
132
133   fAcceptanceCuts = new AliCFAcceptanceCuts("Acceptance", "MC Acceptance Cuts");
134   fAcceptanceCuts->SetMinNHitITS(3);
135   fAcceptanceCuts->SetMinNHitTPC(2);
136   fAcceptanceCuts->SetMinNHitTRD(0);
137
138   // create additional histos for testing purpose
139   fOutput = new AliHFEcollection("histos", "QA histograms");
140   fOutput->CreateTH1F("hNtracks","Number of Tracks per Event", 100, 0, 100);
141   fOutput->CreateTH1F("hPt","Pt of the tracks", 40, 0.1, 10, 0);
142   fOutput->CreateTH1F("kinkIndex", "Kink Index", 400, -200, 200);
143   fOutput->CreateTH1F("itspixel", "ITS PIXEL", 2, 0, 2);
144   fOutput->CreateTH1F("mcmother", "Mother PDG", 1000, -500, 500);
145   fOutput->CreateTH2F("ptres", "Pt Resolution", 40, 0.1, 10, 200, -1.5, 1.5, 0);
146 }
147
148 void AliHFEefficiency::UserExec(Option_t *){
149   //
150   // Event Loop
151   // Filter track, fill Efficiency container
152   //
153   fEfficiency->NewEvent();
154   fFilter->SetRecEvent(fInputEvent);
155   if(fMCEvent){
156     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
157     if(!mcH->InitOk()) return;
158     if(!mcH->TreeK()) return;
159     if(!mcH->TreeTR()) return;
160     fFilter->SetMC(fMCEvent);
161     FilterMC();
162   }
163   fFilter->FilterTracks(dynamic_cast<AliESDEvent *>(fInputEvent));
164   TObjArray *tracks = fFilter->GetFilteredTracks();
165   TIterator *iter = tracks->MakeIterator();
166   fOutput->Fill("hNtracks", tracks->GetEntriesFast());
167   AliESDtrack *track = NULL;
168   while((track = dynamic_cast<AliESDtrack *>(iter->Next()))){
169     fOutput->Fill("hPt", track->Pt());  
170     fOutput->Fill("kinkIndex", track->GetKinkIndex(0));
171     if(TESTBIT(track->GetITSClusterMap(), 0))
172       fOutput->Fill("itspixel",1);
173     else
174       fOutput->Fill("itspixel",0);
175     if(fMCEvent){
176       AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
177       if(mctrack){
178         Int_t motherLabel = mctrack->Particle()->GetFirstMother();
179         if(motherLabel){
180           AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherLabel));
181           if(mother)fOutput->Fill("mcmother", mother->Particle()->GetPdgCode());
182         }
183         fOutput->Fill("ptres", mctrack->Pt(), (track->Pt() - mctrack->Pt())/mctrack->Pt());
184       }
185     }
186   }
187   delete iter;
188   fFilter->Flush();
189   PostData(1, fEfficiency);
190   PostData(2, fOutput->GetList());
191 }
192
193 void AliHFEefficiency::Terminate(Option_t *){
194   //
195   // Evaluate Results
196   //
197   fEfficiency = dynamic_cast<AliHFEcontainer *>(GetOutputData(1));
198   if(!fEfficiency){
199     AliError("No Output data available");
200     return;
201   }
202
203   if(!IsRunningTerminate()) return;
204   PostProcess();
205
206   TList *l = dynamic_cast<TList *>(GetOutputData(2));
207   DrawPtResolution(l);
208 }
209
210 void AliHFEefficiency::FilterMC(){
211   //
212   // Cut MC tracks
213   //
214   Double_t cont[4] = {0., 0., 0., 0.};
215   AliMCParticle *track = NULL;
216   for(Int_t itrack = 0; itrack < fMCEvent->GetNumberOfTracks(); itrack++){
217     track = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrack));
218     if(!fMCcut->IsSelected(track)) continue;
219     cont[0] = track->Pt();
220     cont[1] = track->Eta();
221     cont[2] = track->Phi();
222     cont[3] = track->Charge()/3;
223     //AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(track->Particle()->GetFirstMother()));
224     //if(TMath::Abs(mother->Particle()->GetPdgCode()) != 443) continue;
225     fEfficiency->FillCFContainer("MCFilter", 0, cont);
226     if(!fAcceptanceCuts->IsSelected(track)) continue;
227     fEfficiency->FillCFContainer("MCFilter", 1, cont);
228   }
229 }
230
231 void AliHFEefficiency::Load(const char* filename){
232   //
233   // Load results for post processing
234   //
235   TFile *input = TFile::Open(filename);
236   AliHFEcontainer *cin = dynamic_cast<AliHFEcontainer *>(input->Get("Efficiency"));
237   fEfficiency = dynamic_cast<AliHFEcontainer *>(cin->Clone());
238   input->Close();
239   delete input;
240 }
241
242 void AliHFEefficiency::PostProcess(){
243   //
244   // Calculate Efficiencies
245   //
246   fEfficiency->Print();
247
248   // Prepare containers
249   AliCFContainer *crec  = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signal");
250   AliCFContainer *cmc   = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signalMC");
251   AliCFEffGrid *gridRec = new AliCFEffGrid("effrec", "Efficiency Calculation using rec vars", *crec);
252   AliCFEffGrid *gridMC  = new AliCFEffGrid("effmc", "Efficiency Calculation using mc vars", *cmc);
253
254   TCanvas *ctmp = NULL;
255   for(Int_t ivar = 0; ivar < crec->GetNVar(); ivar++){
256     ctmp = new TCanvas(Form("cEffSignal%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s)", crec->GetVarTitle(ivar)));
257     ctmp->cd();
258     DrawSignalEfficiency(gridRec, crec, ivar);
259     
260     ctmp = new TCanvas(Form("cCutEff%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s)", crec->GetVarTitle(ivar)));
261     ctmp->cd();
262     DrawCutEfficiency(gridRec, crec, ivar);
263   
264     ctmp = new TCanvas(Form("cEffSignalMC%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
265     ctmp->cd();
266     DrawSignalEfficiency(gridMC, cmc, ivar);
267
268     ctmp = new TCanvas(Form("cCutEffMC%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
269     ctmp->cd();
270     DrawCutEfficiency(gridMC, cmc, ivar);
271   } 
272
273   CalculatePTsmearing();
274
275   delete gridRec;
276   delete gridMC;
277   delete crec;
278   delete cmc;
279 }
280
281 void AliHFEefficiency::DrawSignalEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
282   //
283   // Efficiency of a cut step with respect to the signal
284   //
285   TH1 *hEff = NULL;
286   Bool_t kFirst = kTRUE;
287   TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
288   for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
289     eff->CalculateEfficiency(istep, 0);
290     hEff = eff->Project(var);
291     hEff->SetTitle(Form("Signal Efficiency (%s)", cont->GetVarTitle(var)));
292     hEff->GetYaxis()->SetRangeUser(0., 1.6);
293     hEff->GetYaxis()->SetTitle("Efficiency");
294     hEff->SetStats(kFALSE);
295     hEff->SetMarkerStyle(20+istep);
296     hEff->SetMarkerColor(kBlue - 5 + istep);
297     hEff->Draw(kFirst ? "ep" : "epsame");
298     kFirst = kFALSE;
299     leg->AddEntry(hEff, cont->GetStepTitle(istep), "p");
300   }
301   leg->SetBorderSize(0);
302   leg->SetFillStyle(0);
303   leg->Draw();
304   gPad->Update();
305 }
306
307 void AliHFEefficiency::DrawCutEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
308   //
309   // Efficiency of a cut step with respect to the previous one
310   //
311   Bool_t kFirst = kTRUE;
312   TH1 *hEff = NULL;
313   TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
314   TString effTitle;
315   Int_t start = 0, length = 0;
316   for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
317     eff->CalculateEfficiency(istep, istep -1);
318     hEff = eff->Project(var);
319     effTitle = hEff->GetTitle();
320     start = effTitle.Index("projected"); length = effTitle.Length() - start;
321     effTitle.Remove(start - 1, length + 1);
322     effTitle.ReplaceAll("Efficiency: ", "");
323     hEff->SetTitle(Form("Cut Efficiency (%s)", cont->GetVarTitle(var)));
324     hEff->GetYaxis()->SetRangeUser(0., 1.6);
325     hEff->GetYaxis()->SetTitle("Efficiency");
326     hEff->SetStats(kFALSE);
327     hEff->SetMarkerStyle(20+istep);
328     hEff->SetMarkerColor(kBlue - 5 + istep);
329     hEff->Draw(kFirst ? "ep" : "epsame");
330     kFirst = kFALSE;
331     leg->AddEntry(hEff, effTitle.Data(), "p");
332   }
333   leg->SetBorderSize(0);
334   leg->SetFillStyle(0);
335   leg->Draw();
336   gPad->Update();
337 }
338
339 void AliHFEefficiency::CalculatePTsmearing(){
340   //
341   // Efficiency of a cut step with respect to the same cut step filled with MC 
342   // Information
343   //
344   AliCFContainer *cont = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "testfilter_container_signal:testfilter_container_signalMC");
345   AliCFEffGrid *grid = new AliCFEffGrid("effc", "Efficiency Calculation", *cont);
346   Bool_t kFirst = kTRUE;
347   TH1 *hEff = NULL;
348   TCanvas *c2 = new TCanvas("cSmear", "pT modification");
349   c2->cd();
350   TLegend *leg1 = new TLegend(0.5, 0.7, 0.89, 0.89);
351   Int_t nStep = cont->GetNStep()/2;
352   for(Int_t istep = 0; istep < nStep; istep++){
353     grid->CalculateEfficiency(istep, istep + nStep);
354     hEff = grid->Project(0);
355     hEff->SetStats(kFALSE);
356     hEff->SetMarkerStyle(20+istep);
357     hEff->SetMarkerColor(kBlue - 5 + istep);
358     hEff->Draw(kFirst ? "ep" : "epsame");
359     kFirst = kFALSE;
360     leg1->AddEntry(hEff, cont->GetStepTitle(istep), "p");
361   }
362   leg1->Draw();
363   gPad->Update();
364
365   delete cont;
366   delete grid;
367 }
368
369 void AliHFEefficiency::DrawPtResolution(const TList * const l){
370   //
371   // Draw pt resolution
372   //
373   TH2 *h2 = dynamic_cast<TH2 *>(l->FindObject("ptres"));
374   TGraphErrors *mean = new TGraphErrors(h2->GetNbinsX());
375   TGraphErrors *rms = new TGraphErrors(h2->GetNbinsX());
376
377   Double_t pt = 0., pterr = 0., mymean = 0., myrms = 0., mymeanerr = 0., myrmserr = 0.;
378   TH1 *htmp = NULL;
379   for(Int_t imom = 1; imom <= h2->GetXaxis()->GetLast(); imom++){
380     pt = h2->GetXaxis()->GetBinCenter(imom);
381     pterr = h2->GetXaxis()->GetBinWidth(imom)/2.;
382     htmp = h2->ProjectionY("py", imom, imom);
383     mymean = htmp->GetMean();
384     myrms = htmp->GetRMS();
385     mymeanerr = htmp->GetMeanError();
386     myrmserr = htmp->GetRMSError();
387     delete htmp;
388     mean->SetPoint(imom -1, pt, mymean);
389     rms->SetPoint(imom -1, pt, myrms);
390     mean->SetPointError(imom-1,pterr,mymeanerr);
391     rms->SetPointError(imom-1,pterr,myrmserr);
392   }
393
394   TCanvas *c1 = new TCanvas("cPtShift", "pT Shift");
395   c1->cd();
396   mean->SetMarkerStyle(22);
397   mean->SetMarkerColor(kBlue);
398   mean->Draw("ape");
399   rms->SetMarkerStyle(22);
400   rms->SetMarkerColor(kBlack);
401   rms->Draw("pesame");
402   TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
403   leg->AddEntry(mean, "Mean", "lp");
404   leg->AddEntry(rms, "RMS", "lp");
405   leg->Draw();
406   gPad->Update();
407
408