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