Cleanup the code. Fix memory leak. Now inherit from AliAnalysisTaskSE (Antoine, Phili...
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEextraCuts.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 /* $Id$ */
17
18 //
19 // Extra cuts implemented by the ALICE Heavy Flavour Electron Group
20 // Cuts stored here:
21 // - ITS pixels
22 // - TPC cluster ratio
23 // - TRD tracklets
24 //
25 // Authors:
26 //   Markus Fasel <M.Fasel@gsi.de>
27 //
28 #include <TBits.h>
29 #include <TClass.h>
30 #include <TH1F.h>
31 #include <TH2F.h>
32 #include <TList.h>
33 #include <TString.h>
34 #include <TMath.h>
35
36 #include "AliAODTrack.h"
37 #include "AliAODPid.h"
38 #include "AliESDEvent.h"
39 #include "AliESDVertex.h"
40 #include "AliESDtrack.h"
41 #include "AliLog.h"
42 #include "AliMCParticle.h"
43 #include "AliVEvent.h"
44 #include "AliVTrack.h"
45 #include "AliVParticle.h"
46 #include "AliVertexerTracks.h"
47 #include "AliVVertex.h"
48
49 #include "AliHFEextraCuts.h"
50
51 ClassImp(AliHFEextraCuts)
52
53 //______________________________________________________
54 AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
55   AliCFCutBase(name, title),
56   fEvent(NULL),
57   fCutCorrelation(0),
58   fRequirements(0),
59   fMinNClustersTPC(0),
60   fClusterRatioTPC(0.),
61   fMinTrackletsTRD(0),
62   fPixelITS(0),
63   fTOFpid(kFALSE),
64   fTPCclusterDef(0),
65   fTPCclusterRatioDef(0),
66   fCheck(kFALSE),
67   fQAlist(0x0) ,
68   fDebugLevel(0)
69 {
70   //
71   // Default Constructor
72   //
73   memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
74 }
75
76 //______________________________________________________
77 AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
78   AliCFCutBase(c),
79   fEvent(c.fEvent),
80   fCutCorrelation(c.fCutCorrelation),
81   fRequirements(c.fRequirements),
82   fMinNClustersTPC(c.fMinNClustersTPC),
83   fClusterRatioTPC(c.fClusterRatioTPC),
84   fMinTrackletsTRD(c.fMinTrackletsTRD),
85   fPixelITS(c.fPixelITS),
86   fTOFpid(c.fTOFpid),
87   fTPCclusterDef(c.fTPCclusterDef),
88   fTPCclusterRatioDef(c.fTPCclusterRatioDef),
89   fCheck(c.fCheck),
90   fQAlist(0x0),
91   fDebugLevel(0)
92 {
93   //
94   // Copy constructor
95   // Performs a deep copy
96   //
97   memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
98   if(IsQAOn()){
99     fIsQAOn = kTRUE;
100     fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
101     if(fQAlist) fQAlist->SetOwner();
102   }
103 }
104
105 //______________________________________________________
106 AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
107   //
108   // Assignment operator
109   //
110   if(this != &c){
111     AliCFCutBase::operator=(c);
112     fEvent = c.fEvent;
113     fCutCorrelation = c.fCutCorrelation;
114     fRequirements = c.fRequirements;
115     fClusterRatioTPC = c.fClusterRatioTPC;
116     fMinNClustersTPC = c.fMinNClustersTPC;
117     fMinTrackletsTRD = c.fMinTrackletsTRD;
118     fPixelITS = c.fPixelITS;
119     fTPCclusterDef = c.fTPCclusterDef;
120     fTPCclusterRatioDef = c.fTPCclusterRatioDef;
121     fTOFpid = c.fTOFpid;
122     fCheck = c.fCheck;
123     fDebugLevel = c.fDebugLevel;
124     memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
125     if(IsQAOn()){
126       fIsQAOn = kTRUE;
127       fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
128       if(fQAlist) fQAlist->SetOwner();
129     }else fQAlist = 0x0;
130   }
131   return *this;
132 }
133
134 //______________________________________________________
135 AliHFEextraCuts::~AliHFEextraCuts(){
136   //
137   // Destructor
138   //
139   if(fQAlist) delete fQAlist;
140 }
141
142 //______________________________________________________
143 void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
144   //
145   // Set Virtual event an make a copy
146   //
147   if (!event) {
148     AliError("Pointer to AliVEvent !");
149     return;
150   }
151   TString className(event->ClassName());
152   if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
153     AliError("argument must point to an AliESDEvent or AliAODEvent !");
154     return ;
155   }
156   fEvent = (AliVEvent*) event;
157
158 }
159
160 //______________________________________________________
161 Bool_t AliHFEextraCuts::IsSelected(TObject *o){
162   //
163   // Steering function for the track selection
164   //
165   TString type = o->IsA()->GetName();
166   AliDebug(2, Form("Object type %s", type.Data()));
167   if(!type.CompareTo("AliESDtrack") || !type.CompareTo("AliAODTrack")){
168     return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
169   }
170   return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
171 }
172
173 //______________________________________________________
174 Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
175   //
176   // Checks cuts on reconstructed tracks
177   // returns true if track is selected
178   // QA histograms are filled before track selection and for
179   // selected tracks after track selection
180   //
181   AliDebug(1, "Called");
182   ULong64_t survivedCut = 0;    // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
183   if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
184   // Apply cuts
185   Float_t impactR, impactZ;
186   Double_t hfeimpactR, hfeimpactnsigmaR;
187   Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
188   Bool_t tofstep = kTRUE;
189   GetImpactParameters(track, impactR, impactZ);
190   if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
191     // Protection for PbPb
192     GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
193     GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
194   }
195   UInt_t nclsTPC = GetTPCncls(track);
196   // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
197   Double_t ratioTPC = GetTPCclusterRatio(track);
198   UChar_t trdTracklets;
199   trdTracklets = GetTRDnTrackletsPID(track);
200   UChar_t itsPixel = track->GetITSClusterMap();
201   Int_t status1 = GetITSstatus(track, 0);
202   Int_t status2 = GetITSstatus(track, 1);
203   Bool_t statusL0 = CheckITSstatus(status1);
204   Bool_t statusL1 = CheckITSstatus(status2);
205   if(TESTBIT(fRequirements, kMinImpactParamR)){
206     // cut on min. Impact Parameter in Radial direction
207     if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
208   }
209   if(TESTBIT(fRequirements, kMinImpactParamZ)){
210     // cut on min. Impact Parameter in Z direction
211     if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
212   }
213   if(TESTBIT(fRequirements, kMaxImpactParamR)){
214     // cut on max. Impact Parameter in Radial direction
215     if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
216   }
217   if(TESTBIT(fRequirements, kMaxImpactParamZ)){
218     // cut on max. Impact Parameter in Z direction
219     if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
220   }
221   if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
222     // cut on min. HFE Impact Parameter in Radial direction
223     if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
224   }
225   if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
226     // cut on max. HFE Impact Parameter n sigma in Radial direction
227     if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
228   }
229   if(TESTBIT(fRequirements, kClusterRatioTPC)){
230     // cut on min ratio of found TPC clusters vs findable TPC clusters
231     if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
232   }
233   if(TESTBIT(fRequirements, kMinTrackletsTRD)){
234     // cut on minimum number of TRD tracklets
235     AliDebug(1, Form("Min TRD cut: [%d|%d]\n", fMinTrackletsTRD, trdTracklets));
236     if(trdTracklets >= fMinTrackletsTRD) SETBIT(survivedCut, kMinTrackletsTRD);
237   }
238   if(TESTBIT(fRequirements, kMinNClustersTPC)){
239     // cut on minimum number of TRD tracklets
240     AliDebug(1, Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
241     if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
242   }
243   if(TESTBIT(fRequirements, kPixelITS)){
244     // cut on ITS pixel layers
245     AliDebug(1, "ITS cluster Map: ");
246     //PrintBitMap(itsPixel);
247     switch(fPixelITS){
248       case kFirst:
249         AliDebug(2, "First");
250               if(TESTBIT(itsPixel, 0)) 
251                 SETBIT(survivedCut, kPixelITS);
252         else
253                 if(fCheck && !statusL0)
254                   SETBIT(survivedCut, kPixelITS);
255                     break;
256       case kSecond: 
257         AliDebug(2, "Second");
258               if(TESTBIT(itsPixel, 1))
259                 SETBIT(survivedCut, kPixelITS);
260         else
261                 if(fCheck && !statusL1)
262                   SETBIT(survivedCut, kPixelITS);
263                     break;
264       case kBoth: 
265         AliDebug(2, "Both");
266               if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
267                       SETBIT(survivedCut, kPixelITS);
268               else
269           if(fCheck && !(statusL0 && statusL1)) 
270                         SETBIT(survivedCut, kPixelITS);
271               break;
272       case kAny: 
273         AliDebug(2, "Any");
274               if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
275                 SETBIT(survivedCut, kPixelITS);
276         else
277                 if(fCheck && !(statusL0 || statusL1))
278                     SETBIT(survivedCut, kPixelITS);
279                     break;
280       default: 
281         AliDebug(2, "None");
282         break;
283     }
284     AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
285   }
286   if(fTOFpid){
287     // cut on TOF pid
288     tofstep = kFALSE;
289     if(track->GetStatus() & AliESDtrack::kTOFpid) tofstep = kTRUE;
290   
291   }
292   if(fRequirements == survivedCut && tofstep){
293     //
294     // Track selected
295     //
296     AliDebug(2, "Track Survived cuts\n");
297     if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
298     return kTRUE;
299   }
300   AliDebug(2, "Track cut");
301   if(IsQAOn()) FillCutCorrelation(survivedCut);
302   return kFALSE;
303 }
304
305 //______________________________________________________
306 Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
307   //
308   // Checks cuts on Monte Carlo tracks
309   // returns true if track is selected
310   // QA histograms are filled before track selection and for
311   // selected tracks after track selection
312   //
313   return kTRUE; // not yet implemented
314 }
315
316 //______________________________________________________
317 void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
318   //
319   // Fill the QA histograms for ESD tracks
320   // Function can be called before cuts or after cut application (second argument)
321   //
322   const Int_t kNhistos = 6;
323   Float_t impactR, impactZ;
324   GetImpactParameters(track, impactR, impactZ);
325   TH1 *htmp = NULL;
326   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * kNhistos)))) htmp->Fill(impactR);
327   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * kNhistos)))) htmp->Fill(impactZ);
328   // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
329   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * kNhistos)))) htmp->Fill(GetTPCclusterRatio(track));
330   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * kNhistos)))) htmp->Fill(GetTRDnTrackletsPID(track));
331   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * kNhistos)))) htmp->Fill(GetTPCncls(track));
332   UChar_t itsPixel = track->GetITSClusterMap();
333   TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * kNhistos));
334   //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
335   if(pixelHist){
336     Double_t firstEntry = 0.5;
337     if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
338       pixelHist->Fill(firstEntry + 3);
339     else{
340       if(itsPixel & BIT(0)){
341         pixelHist->Fill(firstEntry);
342         if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
343         else pixelHist->Fill(firstEntry + 4);
344       }
345       if(itsPixel & BIT(1)){
346         pixelHist->Fill(firstEntry + 1);
347         if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
348       }
349     }
350   }
351 }
352
353 // //______________________________________________________
354 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
355 //   //
356 //   // Fill the QA histograms for MC tracks
357 //   // Function can be called before cuts or after cut application (second argument)
358 //   // Not yet implemented
359 //   //
360 // }
361
362 //______________________________________________________
363 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
364   //
365   // Fill cut correlation histograms for tracks that didn't pass cuts
366   //
367   const Int_t kNhistos = 6;
368   TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * kNhistos));
369   if(correlation){
370     for(Int_t icut = 0; icut < kNcuts; icut++){
371       if(!TESTBIT(fRequirements, icut)) continue;
372       for(Int_t jcut = icut; jcut < kNcuts; jcut++){
373         if(!TESTBIT(fRequirements, jcut)) continue;
374         if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
375                 correlation->Fill(icut, jcut);
376       }
377     }
378   }
379 }
380
381 //______________________________________________________
382 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
383   //
384   // Add QA histograms
385   // For each cut a histogram before and after track cut is created
386   // Histos before respectively after cut are stored in different lists
387   // Additionally a histogram with the cut correlation is created and stored
388   // in the top directory
389   //
390
391   const Int_t kNhistos = 6;
392   TH1 *histo1D = 0x0;
393   TH2 *histo2D = 0x0;
394   TString cutstr[2] = {"before", "after"};
395
396   if(!fQAlist) fQAlist = new TList;  // for internal representation, not owner
397   for(Int_t icond = 0; icond < 2; icond++){
398     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * kNhistos);
399     fQAlist->AddAt(histo1D, 0 + icond * kNhistos);
400     histo1D->GetXaxis()->SetTitle("Impact Parameter");
401     histo1D->GetYaxis()->SetTitle("Number of Tracks");
402     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * kNhistos);
403     fQAlist->AddAt(histo1D, 1 + icond * kNhistos);
404     histo1D->GetXaxis()->SetTitle("Impact Parameter");
405     histo1D->GetYaxis()->SetTitle("Number of Tracks");
406     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1)), 2 + icond * kNhistos);
407     fQAlist->AddAt(histo1D, 2 + icond * kNhistos);
408     histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
409     histo1D->GetYaxis()->SetTitle("Number of Tracks");
410     qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * kNhistos);
411     fQAlist->AddAt(histo1D, 3 + icond * kNhistos);
412     histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
413     histo1D->GetYaxis()->SetTitle("Number of Tracks");
414     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * kNhistos);
415     fQAlist->AddAt(histo1D, 4 + icond * kNhistos);
416     histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
417     histo1D->GetYaxis()->SetTitle("Number of Tracks");
418     qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * kNhistos);
419     fQAlist->AddAt(histo1D, 5 + icond * kNhistos);
420     histo1D->GetXaxis()->SetTitle("ITS Pixel");
421     histo1D->GetYaxis()->SetTitle("Number of Tracks");
422     Int_t first = histo1D->GetXaxis()->GetFirst();
423     TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
424     for(Int_t ilabel = 0; ilabel < 6; ilabel++)
425       histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
426   }
427   // Add cut correlation
428   qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * kNhistos);
429   fQAlist->AddAt(histo2D, 2 * kNhistos);
430   TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel", "kMinHFEImpactParamR", "kMinHFEImpactParamNsigmaR", "TPC Number of clusters"};
431   Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
432   for(Int_t icut = 0; icut < kNcuts; icut++){
433     histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
434     histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
435   }
436 }
437
438 //______________________________________________________
439 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
440   for(Int_t ibit = 32; ibit--; )
441     printf("%d", bitmap & BIT(ibit) ? 1 : 0);
442   printf("\n");
443 }
444
445 //______________________________________________________
446 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
447   //
448   // Check whether ITS area is dead
449   //
450   Bool_t status;
451   switch(itsStatus){
452     case 2: status = kFALSE; break;
453     case 3: status = kFALSE; break;
454     case 7: status = kFALSE; break;
455     default: status = kTRUE;
456   }
457   return status;
458 }
459
460 //______________________________________________________
461 Int_t AliHFEextraCuts::GetTRDnTrackletsPID(AliVTrack *track){
462         //
463         // Get Number of TRD tracklets
464         //
465         Int_t nTracklets = 0;
466         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
467                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
468                 if(esdtrack) nTracklets = esdtrack->GetTRDntrackletsPID();
469         } else if(!TString(track->IsA()->GetName()).CompareTo("AliAODTrack")){
470                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
471                 AliAODPid *pidobject = NULL;
472     if(aodtrack) pidobject = aodtrack->GetDetPid();
473                 // this is normally NOT the way to do this, but due to limitation in the
474                 // AOD track it is not possible in a different way
475                 if(pidobject){
476                         Float_t *trdmom = pidobject->GetTRDmomentum();
477                         for(Int_t ily = 0; ily < 6; ily++){
478                                 if(trdmom[ily] > -1) nTracklets++;
479                         }
480                 } else nTracklets = 6;  // No Cut possible
481         }
482         return nTracklets;
483 }
484
485 //______________________________________________________
486 Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
487         //
488         // Check ITS layer status
489         //
490         Int_t status = 0;
491         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
492                 Int_t det;
493                 Float_t xloc, zloc;
494                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
495                 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
496         }
497         return status;
498 }
499
500 //______________________________________________________
501 UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
502         //
503         // Get Number of findable clusters in the TPC
504         //
505         Int_t nClusters = 0; // in case no Information available consider all clusters findable
506         TString type = track->IsA()->GetName();
507         if(!type.CompareTo("AliESDtrack")){
508     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
509     if(esdtrack){ // coverity
510       if(TESTBIT(fTPCclusterDef, kFoundIter1)){
511                     nClusters = esdtrack->GetTPCNclsIter1();
512         AliDebug(2, ("Using def kFoundIter1"));
513       } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
514         AliDebug(2, ("Using def kCrossedRows"));
515         nClusters = static_cast<UInt_t>(esdtrack->GetTPCClusterInfo(2,1));
516       } else{
517         AliDebug(2, ("Using def kFound"));
518                     nClusters = esdtrack->GetTPCNcls();
519       }
520     }
521   }
522         else if(!type.CompareTo("AliAODTrack")){
523                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
524     if(aodtrack){
525             const TBits &tpcmap = aodtrack->GetTPCClusterMap();
526             for(UInt_t ibit = 0; ibit < tpcmap.GetNbits(); ibit++)
527                     if(tpcmap.TestBitNumber(ibit)) nClusters++;
528     }
529         }
530         return nClusters;
531 }
532
533 //______________________________________________________
534 Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
535   //
536   // Get Ratio of found / findable clusters for different definitions
537   // Only implemented for ESD tracks
538   //
539   Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
540   TString type = track->IsA()->GetName();
541   if(!type.CompareTo("AliESDtrack")){
542     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
543     if(esdtrack){ // coverity
544       if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
545         AliDebug(2, "Using ratio def kCROverFindable");
546         clusterRatio = esdtrack->GetTPCNclsF() ? esdtrack->GetTPCClusterInfo(2,1)/static_cast<Double_t>(esdtrack->GetTPCNclsF()) : 1.; // crossed rows/findable
547       } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
548         AliDebug(2, "Using ratio def kFoundOverCR");
549         clusterRatio = esdtrack->GetTPCClusterInfo(2,0); // found/crossed rows
550       } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
551         AliDebug(2, "Using ratio def kFoundOverFindableIter1");
552         clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
553       } else {
554         AliDebug(2, "Using ratio def kFoundOverFindable");
555         clusterRatio = esdtrack->GetTPCNclsF() ? static_cast<Double_t>(esdtrack->GetTPCNcls())/static_cast<Double_t>(esdtrack->GetTPCNclsF()) : 1.; // found/findable
556       }
557     }
558   }
559   return clusterRatio;
560 }
561
562 //______________________________________________________
563 void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
564         //
565         // Get impact parameter
566         //
567         TString type = track->IsA()->GetName();
568         if(!type.CompareTo("AliESDtrack")){
569                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
570     if(esdtrack) esdtrack->GetImpactParameters(radial, z);
571         }
572         else if(!type.CompareTo("AliAODTrack")){
573                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
574     if(aodtrack){
575                   Double_t xyz[3];
576                   aodtrack->XYZAtDCA(xyz);
577                   z = xyz[2];
578                   radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
579     }
580         }
581 }
582
583 //______________________________________________________
584 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
585         //
586         // Get HFE impact parameter (with recalculated primary vertex)
587         //
588         dcaxy=0;
589         dcansigmaxy=0;
590   if(!fEvent){
591     AliDebug(1, "No Input event available\n");
592     return;
593   }
594   const Double_t kBeampiperadius=3.;
595   TString type = track->IsA()->GetName();
596   Double_t dca[2]={-999.,-999.};
597   Double_t cov[3]={-999.,-999.,-999.};
598
599   // recalculate primary vertex
600   AliVertexerTracks vertexer(fEvent->GetMagneticField());
601   vertexer.SetITSMode();
602   vertexer.SetMinClusters(4);
603         Int_t skipped[2];
604   skipped[0] = track->GetID();
605   vertexer.SetSkipTracks(1,skipped);
606   AliVVertex *vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
607   vertexer.SetSkipTracks(1,skipped);
608   if(vtxESDSkip->GetNContributors()<2) return;
609
610   // Getting the DCA
611   // Propagation always done on a working copy to not disturb the track params of the original track
612   AliESDtrack *esdtrack = NULL;
613   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
614     // Case ESD track: take copy constructor
615     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
616     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
617   } else {
618     // Case AOD track: take different constructor
619     esdtrack = new AliESDtrack(track);
620   }
621   if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dca, cov)){
622     // protection
623     dcaxy = dca[0];
624     if(cov[0]) dcansigmaxy = dcaxy/TMath::Sqrt(cov[0]);
625     if(!cov[0]) dcansigmaxy = -99.;
626   }
627   delete esdtrack;
628   delete vtxESDSkip;
629 }
630
631
632 //______________________________________________________
633 void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
634         //
635         // Get HFE impact parameter cut(pt dependent)
636         //
637   
638   TString type = track->IsA()->GetName();
639   if(!type.CompareTo("AliESDtrack")){
640     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
641     if(!esdtrack) return;
642     Double_t pt = esdtrack->Pt();       
643     //hfeimpactRcut=0.0064+0.078*exp(-0.56*pt);  // used Carlo's old parameter 
644     hfeimpactRcut=0.011+0.077*exp(-0.65*pt); // used Carlo's new parameter
645     hfeimpactnsigmaRcut=3; // 3 sigma trail cut
646   }
647 }
648