Cleanup the code. Fix memory leak. Now inherit from AliAnalysisTaskSE (Antoine, Phili...
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEtrackFilter.cxx
CommitLineData
70da6c5a 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**************************************************************************/
27de2dfb 15
16/* $Id$ */
17
70da6c5a 18//
19// Track filter class
20// Apply cut steps to all tracks in one event and returns a list of
21// filtered tracks
22//
23// Author:
24// Markus Fasel <M.Fasel@gsi.de>
25//
26#include <TMath.h>
27#include <TString.h>
28
29#include "AliAnalysisCuts.h"
30#include "AliCFContainer.h"
31#include "AliCFParticleGenCuts.h"
32#include "AliCFTrackKineCuts.h"
33#include "AliCFTrackIsPrimaryCuts.h"
34#include "AliCFTrackQualityCuts.h"
35#include "AliESDEvent.h"
36#include "AliESDtrack.h"
37#include "AliMCEvent.h"
38#include "AliMCParticle.h"
39
40#include "AliHFEcontainer.h"
41#include "AliHFEcutStep.h"
42#include "AliHFEextraCuts.h"
43#include "AliHFEtrackFilter.h"
44#include "AliHFEtools.h"
45
46ClassImp(AliHFEtrackFilter)
47
48//__________________________________________________________________
49AliHFEtrackFilter::AliHFEtrackFilter(const Char_t *name) :
50 TNamed(name, ""),
51 fFilteredTracks(NULL),
52 fCutSteps(NULL),
53 fEfficiencyContainers(NULL),
54 fMC(NULL),
55 fMCsignal(NULL),
56 fPtBins(0),
57 fEtaBins(0),
58 fPhiBins(0),
59 fPtBinning(NULL),
60 fEtaBinning(NULL),
61 fPhiBinning(NULL)
62{
63 //
64 // Constructor
65 //
66 fCutSteps = new TObjArray;
67 fFilteredTracks = new TObjArray;
68 fEfficiencyContainers = new TObjArray(4);
69}
70
71//__________________________________________________________________
72AliHFEtrackFilter::AliHFEtrackFilter(const AliHFEtrackFilter &ref):
73 TNamed(ref),
74 fFilteredTracks(NULL),
75 fCutSteps(NULL),
76 fEfficiencyContainers(NULL),
77 fMC(NULL),
78 fMCsignal(NULL),
79 fPtBins(0),
80 fEtaBins(0),
81 fPhiBins(0),
82 fPtBinning(NULL),
83 fEtaBinning(NULL),
84 fPhiBinning(NULL)
85{
86 //
87 // Copy constructor
88 //
89 ref.Copy(*this);
90}
91
92//__________________________________________________________________
93AliHFEtrackFilter &AliHFEtrackFilter::operator=(const AliHFEtrackFilter &o){
94 //
95 // Assignment operator
96 //
97 if(this != &o){
98 o.Copy(*this);
99 }
100 return *this;
101}
102
103//__________________________________________________________________
104AliHFEtrackFilter::~AliHFEtrackFilter(){
105 //
106 // Destructor
107 //
108
109 // remove correction framework content
110 if(fPtBinning) delete fPtBinning;
111 if(fEtaBinning) delete fEtaBinning;
112 if(fPhiBinning) delete fPhiBinning;
113 if(TestBit(kOwnCFContainers))
114 fEfficiencyContainers->Delete();
115 delete fEfficiencyContainers;
116 // remove default content
117 delete fFilteredTracks;
118 delete fCutSteps;
119}
120
121//__________________________________________________________________
122void AliHFEtrackFilter::Copy(TObject &o) const{
123 //
124 // Copy content into the object o
125 //
126 TNamed::Copy(o);
127 AliHFEtrackFilter &target = dynamic_cast<AliHFEtrackFilter &>(o);
128
129 // Make copy
130 target.fFilteredTracks = dynamic_cast<TObjArray *>(fFilteredTracks->Clone());
131 target.fCutSteps = dynamic_cast<TObjArray *>(fCutSteps->Clone());
132 target.fEfficiencyContainers = dynamic_cast<TObjArray *>(fEfficiencyContainers->Clone());
133 target.fMC = fMC;
134 target.fMCsignal = dynamic_cast<AliHFEcutStep *>(fMCsignal->Clone());
135
136 target.fPtBins = fPtBins;
137 target.fEtaBins = fEtaBins;
138 target.fPhiBins = fPhiBins;
139 target.fPtBinning = new Double_t[fPtBins];
140 memcpy(target.fPtBinning, fPtBinning, sizeof(Double_t) * fPtBins);
141 target.fEtaBinning = new Double_t[fEtaBins];
142 memcpy(target.fEtaBinning, fEtaBinning, sizeof(Double_t) * fEtaBins);
143 target.fPhiBinning = new Double_t[fPhiBins];
144 memcpy(target.fPhiBinning, fPhiBinning, sizeof(Double_t) * fPhiBins);
145}
146
147//__________________________________________________________________
148void AliHFEtrackFilter::InitCF(){
149 //
150 // Initialize correction framework container
151 // No HFE container given
152 // Only for testing purpose
153 //
154 const Char_t *cnames[4] = {"_container", "_container_signal", "_container_MC", "_container_signalMC"};
155 const Char_t *ctitlesAppends[4] = {"(All Rec)", "(All Signals)", "(All Rec (MC))", "(All Signals(MC))"};
156
157 // Create the binning if not done from outside
158 if(!fPtBins){
159 fPtBins = 40;
160 fPtBinning = AliHFEtools::MakeLogarithmicBinning(fPtBins, 0.1, 10);
161 }
162 if(!fEtaBins){
163 fEtaBins = 8;
164 fEtaBinning = AliHFEtools::MakeLinearBinning(fEtaBins, -0.9, 0.9);
165 }
166 if(!fPhiBins){
167 fPhiBins = 18;
168 fPhiBinning = AliHFEtools::MakeLinearBinning(fPhiBins, 0, 2*TMath::Pi());
169 }
170 Double_t chargeBins[3] = {-1.1, 0., 1.1};
171 Int_t nStep = fCutSteps->GetEntriesFast()+1;
172 Int_t nBins[4] = {fPtBins, fEtaBins, fPhiBins + 2};
173 AliCFContainer *ctmp = NULL;
174 for(Int_t icont = 0; icont < 4; icont++){
175 TString containername = GetName() + TString(cnames[icont]);
176 TString containertitle = TString("Container for filter ") + GetName() + TString(ctitlesAppends[icont]);
177 ctmp = new AliCFContainer(containername.Data(), containertitle.Data(), nStep, 4, nBins);
178 SetBit(kOwnCFContainers, kTRUE);
179
180 // Set the binning
181 ctmp->SetBinLimits(0, fPtBinning);
182 ctmp->SetBinLimits(1, fEtaBinning);
183 ctmp->SetBinLimits(2, fPhiBinning);
184 ctmp->SetBinLimits(3, chargeBins);
185
186 // Label variable names
187 ctmp->SetVarTitle(0, "pt");
188 ctmp->SetVarTitle(1, "eta");
189 ctmp->SetVarTitle(2, "phi");
190 ctmp->SetVarTitle(3, "charge");
191
192 // Label step name
193 ctmp->SetStepTitle(0, "No Cuts");
194 AliHFEcutStep *cutStep = NULL;
195 for(Int_t istep = 0; istep < fCutSteps->GetEntriesFast(); istep++){
196 cutStep = dynamic_cast<AliHFEcutStep *>(fCutSteps->UncheckedAt(istep));
bf892a6a 197 if(cutStep) ctmp->SetStepTitle(istep + 1, cutStep->GetName());
70da6c5a 198 }
199 fEfficiencyContainers->AddAt(ctmp, 0);
200 }
201 OwnContainers();
202}
203
204//__________________________________________________________________
205void AliHFEtrackFilter::InitCF(AliHFEcontainer *cont){
206 //
207 // Initialize Correction Framework container
208 // Use the given HFE container
209 // HFE standard way to use the correction Framework
210 //
211
212 const Char_t *cnames[4] = {"_container", "_container_signal", "_container_MC", "_container_signalMC"};
213 const Char_t *ctitlesAppends[4] = {"(All Rec)", "(All Signals)", "(All Rec (MC))", "(All Signals(MC))"};
214 Int_t nStep = fCutSteps->GetEntriesFast()+1;
215 AliCFContainer *ctmp = NULL;
216 for(Int_t icont = 0; icont < 4; icont++){
217 TString contname = GetName() + TString(cnames[icont]);
218 TString contTitle = TString("Container for filter ") + GetName() + TString(ctitlesAppends[icont]);
219 //printf("Adding container %s: %s\n", contname.Data(), contTitle.Data());
220 cont->CreateContainer(contname.Data(), contTitle.Data(), nStep);
221 fEfficiencyContainers->AddAt((ctmp = cont->GetCFContainer(contname.Data())), icont);
222
223 // Label step name
224 ctmp->SetStepTitle(0, "No Cuts");
225 AliHFEcutStep *cutStep = NULL;
226 for(Int_t istep = 0; istep < fCutSteps->GetEntriesFast(); istep++){
227 cutStep = dynamic_cast<AliHFEcutStep *>(fCutSteps->UncheckedAt(istep));
bf892a6a 228 if(cutStep) ctmp->SetStepTitle(istep + 1, cutStep->GetName());
70da6c5a 229 }
230 }
231 ReleaseContainers();
232}
233
234//__________________________________________________________________
235void AliHFEtrackFilter::FilterTracks(AliESDEvent * const event){
236 //
237 // Perform track filtering
238 // Check each cut step one by one and select tracks which pass
239 // all cuts. If the correction framework is initialized, a correction#
240 // framework container will be filled for each cut step.
241 //
242 AliHFEcutStep *cutStep = NULL;
243 AliESDtrack *track = NULL;
244 AliMCParticle *mctrack = NULL;
245 Double_t cont[4] = {0., 0., 0., 0.}, contMC[4] = {0., 0., 0., 0.};
246 Bool_t goodTrack = kTRUE, signal = kFALSE;
247 Int_t nStep = fCutSteps->GetEntriesFast();
248 AliCFContainer *call = dynamic_cast<AliCFContainer *>(fEfficiencyContainers->At(0));
249 AliCFContainer *cSignal = dynamic_cast<AliCFContainer *>(fEfficiencyContainers->At(1));
250 AliCFContainer *callMC = dynamic_cast<AliCFContainer *>(fEfficiencyContainers->At(2));
251 AliCFContainer *cSignalMC = dynamic_cast<AliCFContainer *>(fEfficiencyContainers->At(3));
252 for(Int_t itrack = 0; itrack < event->GetNumberOfTracks(); itrack++){
253 signal = kFALSE;
254 track = event->GetTrack(itrack);
255 // check Monte-Carlo Information if available
256 if(fMC){
257 mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(track->GetLabel())));
258 if(mctrack){
259 //AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMC->GetTrack(mctrack->Particle()->GetFirstMother()));
260 //AliInfo(Form("Label %d, Mother %d", track->GetLabel(), mother->Particle()->GetPdgCode()));
261 contMC[0] = mctrack->Pt();
262 contMC[1] = mctrack->Eta();
263 contMC[2] = mctrack->Phi();
264 contMC[3] = mctrack->Charge();
265 if(fMCsignal->IsSelected(mctrack)) signal = kTRUE;
266 //if(TMath::Abs(mother->Particle()->GetPdgCode()) != 443) signal = kFALSE;
267 //AliInfo(Form("Signal? %s", signal ? "Yes": "No"));
268 }
269 }
270
271 // Fill Array without cut
272 cont[0] = track->Pt();
273 cont[1] = track->Eta();
274 cont[2] = track->Phi();
275 cont[3] = track->Charge();
276 if(call) call->Fill(cont, 0);
277 if(callMC) callMC->Fill(contMC, 0);
278 if(signal){
279 if(cSignal) cSignal->Fill(cont, 0);
280 if(cSignalMC) cSignalMC->Fill(contMC, 0);
281 }
282 // cut the track
283 goodTrack = kTRUE;
284 for(Int_t icut = 0; icut < nStep; icut++){
285 cutStep = dynamic_cast<AliHFEcutStep *>(fCutSteps->UncheckedAt(icut));
e3ae862b 286 if(cutStep && (!cutStep->IsSelected(track))){
70da6c5a 287 // track cut away
288 goodTrack = kFALSE;
289 break;
290 }
291
292 // Track survived cut step: Fill container
293 if(call) call->Fill(cont, icut + 1);
294 if(callMC) callMC->Fill(contMC, icut + 1);
295 if(signal){
296 if(cSignal) cSignal->Fill(cont, icut + 1);
297 if(cSignalMC) cSignalMC->Fill(contMC, icut + 1);
298 }
299 }
300
301 // Append track to the list of filtered tracks
302 if(goodTrack) fFilteredTracks->Add(track);
303 }
304}
305
306//__________________________________________________________________
307void AliHFEtrackFilter::GenerateCutSteps(){
308 //
309 // Make the default cuts
310 //
311 MakeCutStepRecKineITSTPC();
312 MakeCutStepPrimary();
313 MakeCutStepHFEITS();
314 MakeCutStepHFETRD();
315
316 MakeMCSignalCuts();
317}
318
319//__________________________________________________________________
320void AliHFEtrackFilter::AddCutStep(AliHFEcutStep *step){
321 //
322 // Add new cut step to the filter
323 //
324 if(!fCutSteps->FindObject(step->GetName())) fCutSteps->Add(step);
325}
326
327//__________________________________________________________________
328AliHFEcutStep *AliHFEtrackFilter::GetCutStep(Int_t istep){
329 //
330 // Getter for single cut step
331 //
332 return dynamic_cast<AliHFEcutStep *>(fCutSteps->At(istep));
333}
334
335//__________________________________________________________________
336AliHFEcutStep *AliHFEtrackFilter::GetCutStep(const Char_t *name){
337 //
338 // Getter for single cut step (by step name)
339 //
340 return dynamic_cast<AliHFEcutStep *>(fCutSteps->FindObject(name));
341}
342
343//__________________________________________________________________
344void AliHFEtrackFilter::Flush(){
345 //
346 // Empty track container
347 //
348 fFilteredTracks->Clear();
349}
350
351//__________________________________________________________________
352void AliHFEtrackFilter::SetPtBins(Int_t nBins, Double_t *binning){
353 //
354 // User defined pt binning
355 //
356 fPtBins = nBins;
357 fPtBinning = new Double_t[fPtBins + 1];
358 memcpy(fPtBinning, binning, sizeof(Double_t) * nBins);
359}
360
361//__________________________________________________________________
362void AliHFEtrackFilter::SetEtaBins(Int_t nBins, Double_t *binning){
363 //
364 // User defined eta binning
365 //
366 fEtaBins = nBins;
367 fEtaBinning = new Double_t[fEtaBins + 1];
368 memcpy(fEtaBinning, binning, sizeof(Double_t) * nBins);
369}
370
371//__________________________________________________________________
372void AliHFEtrackFilter::SetPhiBins(Int_t nBins, Double_t *binning){
373 //
374 // User defined phi binning
375 //
376 fPhiBins = nBins;
377 fPhiBinning = new Double_t[fPhiBins + 1];
378 memcpy(fPhiBinning, binning, sizeof(Double_t) * nBins);
379}
380
381//__________________________________________________________________
382AliHFEcutStep *AliHFEtrackFilter::MakeCutStepRecKineITSTPC(){
383 //
384 // Make the cut step for Rec Kine
385 // Cut step is already included in the filter
386 //
387 AliHFEcutStep *fCutStep = new AliHFEcutStep("RecKineITSTPC");
388
faee3b18 389 AliCFTrackQualityCuts *trackQuality = new AliCFTrackQualityCuts((Char_t *)"QualityRec", (Char_t *)"REC Track Quality Cuts");
390 trackQuality->SetMinNClusterTPC(80);
70da6c5a 391 trackQuality->SetMaxChi2PerClusterTPC(3.5);
392 trackQuality->SetStatus(AliESDtrack::kTPCrefit | AliESDtrack::kITSrefit);
393 trackQuality->SetMaxCovDiagonalElements(2., 2., 0.5, 0.5, 2);
394 fCutStep->AddCut(trackQuality);
395
396 AliHFEextraCuts *hfecuts = new AliHFEextraCuts("HFETPC","Extra cuts from the HFE group");
e3ae862b 397 hfecuts->SetClusterRatioTPC(0.6, AliHFEextraCuts::kFoundOverCR);
70da6c5a 398 fCutStep->AddCut(hfecuts);
399
faee3b18 400 AliCFTrackKineCuts *kineCuts = new AliCFTrackKineCuts((Char_t *)"RecKine", (Char_t *)"REC Kine Cuts");
401 kineCuts->SetPtRange(0.1, 20);
402 kineCuts->SetEtaRange(-0.8, 0.8);
70da6c5a 403 fCutStep->AddCut(kineCuts);
404
405 AddCutStep(fCutStep);
406 return fCutStep;
407}
408
409//__________________________________________________________________
410AliHFEcutStep *AliHFEtrackFilter::MakeCutStepPrimary(){
411 //
412 // Make cut on primaries
413 // Cut step is already included in the filter
414 //
415 AliHFEcutStep *fCutStep = new AliHFEcutStep("Primary");
416
faee3b18 417 AliCFTrackIsPrimaryCuts *primaryCut = new AliCFTrackIsPrimaryCuts((Char_t *)"PrimaryCuts", (Char_t *)"REC Primary Cuts");
418 primaryCut->SetMaxDCAToVertexXY(0.5);
419 primaryCut->SetMaxDCAToVertexZ(5.);
70da6c5a 420 primaryCut->SetAcceptKinkDaughters(kFALSE);
421 fCutStep->AddCut(primaryCut);
422
423 AddCutStep(fCutStep);
424 return fCutStep;
425}
426
427//__________________________________________________________________
428AliHFEcutStep *AliHFEtrackFilter::MakeCutStepHFEITS(){
429 //
430 // Add special ITS cuts
431 // Cut step is already included in the filter
432 //
433 AliHFEcutStep *fCutStep = new AliHFEcutStep("HFEITS");
434
faee3b18 435 AliHFEextraCuts *hfecuts = new AliHFEextraCuts((Char_t *)"HFEPixelsCuts",(Char_t *)"Extra cuts from the HFE group");
70da6c5a 436 hfecuts->SetRequireITSpixel(AliHFEextraCuts::kFirst);
437 //hfecuts->SetCheckITSstatus(kTRUE);
438 fCutStep->AddCut(hfecuts);
439
440 AddCutStep(fCutStep);
441 return fCutStep;
442}
443
444//__________________________________________________________________
445AliHFEcutStep *AliHFEtrackFilter::MakeCutStepHFETRD(){
446 //
447 // Add special TRD cut
448 // Cut step is already included in the filter
449 //
450 AliHFEcutStep *fCutStep = new AliHFEcutStep("HFETRD");
451
452 AliHFEextraCuts *hfecuts = new AliHFEextraCuts("HFETRDCuts","Extra cuts from the HFE group");
453 hfecuts->SetMinTrackletsTRD(0);
454 fCutStep->AddCut(hfecuts);
455
456 AddCutStep(fCutStep);
457 return fCutStep;
458}
459
460//__________________________________________________________________
461AliHFEcutStep *AliHFEtrackFilter::MakeMCSignalCuts(){
462 //
463 // Define MC Signal
464 // Cut step is already included in the filter
465 //
466 fMCsignal = new AliHFEcutStep("MCSignal");
faee3b18 467 AliCFParticleGenCuts *genCuts = new AliCFParticleGenCuts((Char_t *)"fCutsGenMC", (Char_t *)"Particle Generation Cuts");
70da6c5a 468 genCuts->SetRequireIsCharged();
469 genCuts->SetRequireIsPrimary();
faee3b18 470 genCuts->SetProdVtxRange2D();
471 genCuts->SetProdVtxRangeX(0, 1);
472 genCuts->SetProdVtxRangeY(0, 1);
70da6c5a 473 genCuts->SetRequirePdgCode(11, kTRUE);
474 fMCsignal->AddCut(genCuts);
475
faee3b18 476 AliCFTrackKineCuts *kineMCcuts = new AliCFTrackKineCuts((Char_t *)"fCutsKineMC",(Char_t *)"MC Kine Cuts");
477 kineMCcuts->SetPtRange(0.1, 20.);
478 kineMCcuts->SetEtaRange(-0.8, 0.8);
70da6c5a 479 fMCsignal->AddCut(kineMCcuts);
480
481 return fMCsignal;
482}
483
484//__________________________________________________________________
485void AliHFEtrackFilter::SetMC(AliMCEvent * const mc){
486 //
487 // Publish MC event to the single cut steps
488 //
489 fMC = mc;
490 fMCsignal->SetMC(fMC);
bf892a6a 491 AliHFEcutStep *cs = NULL;
70da6c5a 492 for(Int_t icut = 0; icut < fCutSteps->GetEntriesFast(); icut++)
bf892a6a 493 if((cs = dynamic_cast<AliHFEcutStep *>(fCutSteps->UncheckedAt(icut)))) cs->SetMC(fMC);
70da6c5a 494}
495
496//__________________________________________________________________
497void AliHFEtrackFilter::SetRecEvent(AliVEvent *rec){
498 //
499 // Publish MC event to the single cut steps
500 //
bf892a6a 501 AliHFEcutStep *cs = NULL;
70da6c5a 502 for(Int_t icut = 0; icut < fCutSteps->GetEntriesFast(); icut++)
bf892a6a 503 if((cs = dynamic_cast<AliHFEcutStep *>(fCutSteps->UncheckedAt(icut)))) cs->SetRecEvent(rec);
70da6c5a 504}
505
506
507//__________________________________________________________________
508AliCFContainer *AliHFEtrackFilter::GetEfficiencyContainer(Int_t icont){
509 //
510 // return EfficiencyContainer
511 //
512 if(icont >= 4 || icont < 0) return NULL;
513 return dynamic_cast<AliCFContainer *>(fEfficiencyContainers->At(icont));
514}