664a5ee02458fe9d5f4588ed37c80918ca30afb3
[u/mrichter/AliRoot.git] / PWGHF / 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 // Extra cuts implemented by the ALICE Heavy Flavour Electron Group
17 // Cuts stored here:
18 // - ITS pixels
19 // - TPC cluster ratio
20 // - TRD tracklets
21 //
22 // Authors:
23 //   Markus Fasel <M.Fasel@gsi.de>
24 //
25 #include <TBits.h>
26 #include <TClass.h>
27 #include <TH1F.h>
28 #include <TH2F.h>
29 #include <TList.h>
30 #include <TString.h>
31 #include <TMath.h>
32
33 #include "AliAODTrack.h"
34 #include "AliAODPid.h"
35 #include "AliAODVertex.h"
36 #include "AliAODTrack.h"
37 #include "AliESDEvent.h"
38 #include "AliESDVertex.h"
39 #include "AliESDtrack.h"
40 #include "AliLog.h"
41 #include "AliMCParticle.h"
42 #include "AliVEvent.h"
43 #include "AliVTrack.h"
44 #include "AliVParticle.h"
45 #include "AliVertexerTracks.h"
46 #include "AliVVertex.h"
47
48 #include "AliHFEextraCuts.h"
49
50 ClassImp(AliHFEextraCuts)
51
52 const Int_t AliHFEextraCuts::fgkNQAhistos = 8;
53
54 //______________________________________________________
55 AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
56   AliCFCutBase(name, title),
57   fEvent(NULL),
58   fCutCorrelation(0),
59   fRequirements(0),
60   fMinNClustersTPC(0),
61   fMinNClustersTPCPID(0),
62   fClusterRatioTPC(0.),
63   fMinTrackletsTRD(0),
64   fMinNbITScls(0),
65   fTRDtrackletsExact(0),
66   fPixelITS(0),
67   fDriftITS(0),
68   fTPCclusterDef(0),
69   fTPCclusterRatioDef(0),
70   fFractionTPCShared(-1.0),
71   fAbsHFEImpactParamNsigmaR(kTRUE),
72   fTOFsignalDx(1.0),
73   fTOFsignalDz(1.0),
74   fCheck(kFALSE),
75   fQAlist(0x0) ,
76   fDebugLevel(0)
77 {
78   //
79   // Default Constructor
80   //
81   //printf("Set the number of min ITS clusters %d\n",(Int_t)fMinNbITScls);
82   memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
83   memset(fIPcutParam, 0, sizeof(Float_t) * 4);
84 }
85
86 //______________________________________________________
87 AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
88   AliCFCutBase(c),
89   fEvent(c.fEvent),
90   fCutCorrelation(c.fCutCorrelation),
91   fRequirements(c.fRequirements),
92   fMinNClustersTPC(c.fMinNClustersTPC),
93   fMinNClustersTPCPID(c.fMinNClustersTPCPID),
94   fClusterRatioTPC(c.fClusterRatioTPC),
95   fMinTrackletsTRD(c.fMinTrackletsTRD),
96   fMinNbITScls(c.fMinNbITScls),
97   fTRDtrackletsExact(c.fTRDtrackletsExact),
98   fPixelITS(c.fPixelITS),
99   fDriftITS(c.fDriftITS),
100   fTPCclusterDef(c.fTPCclusterDef),
101   fTPCclusterRatioDef(c.fTPCclusterRatioDef),
102   fFractionTPCShared(c.fFractionTPCShared),
103   fAbsHFEImpactParamNsigmaR(c.fAbsHFEImpactParamNsigmaR),
104   fTOFsignalDx(c.fTOFsignalDx),
105   fTOFsignalDz(c.fTOFsignalDz),
106   fCheck(c.fCheck),
107   fQAlist(0x0),
108   fDebugLevel(0)
109 {
110   //
111   // Copy constructor
112   // Performs a deep copy
113   //
114   memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
115   memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
116   if(IsQAOn()){
117     fIsQAOn = kTRUE;
118     fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
119     if(fQAlist) fQAlist->SetOwner();
120   }
121 }
122
123 //______________________________________________________
124 AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
125   //
126   // Assignment operator
127   //
128   if(this != &c){
129     AliCFCutBase::operator=(c);
130     fEvent = c.fEvent;
131     fCutCorrelation = c.fCutCorrelation;
132     fRequirements = c.fRequirements;
133     fClusterRatioTPC = c.fClusterRatioTPC;
134     fMinNClustersTPC = c.fMinNClustersTPC;
135     fMinNClustersTPCPID = c.fMinNClustersTPCPID;
136     fMinTrackletsTRD = c.fMinTrackletsTRD;
137     fMinNbITScls = c.fMinNbITScls;
138     fTRDtrackletsExact = c.fTRDtrackletsExact;
139     fPixelITS = c.fPixelITS;
140     fDriftITS = c.fDriftITS;
141     fTPCclusterDef = c.fTPCclusterDef;
142     fTPCclusterRatioDef = c.fTPCclusterRatioDef;
143     fFractionTPCShared = c.fFractionTPCShared;
144     fAbsHFEImpactParamNsigmaR = c.fAbsHFEImpactParamNsigmaR;
145     fTOFsignalDx = c.fTOFsignalDx;
146     fTOFsignalDz = c.fTOFsignalDz;
147     fCheck = c.fCheck;
148     fDebugLevel = c.fDebugLevel;
149     memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
150     memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
151     if(IsQAOn()){
152       fIsQAOn = kTRUE;
153       fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
154       if(fQAlist) fQAlist->SetOwner();
155     }else fQAlist = 0x0;
156   }
157   return *this;
158 }
159
160 //______________________________________________________
161 AliHFEextraCuts::~AliHFEextraCuts(){
162   //
163   // Destructor
164   //
165   if(fQAlist) delete fQAlist;
166 }
167
168 //______________________________________________________
169 void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
170   //
171   // Set Virtual event an make a copy
172   //
173   if (!event) {
174     AliError("Pointer to AliVEvent !");
175     return;
176   }
177   TString className(event->ClassName());
178   if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
179     AliError("argument must point to an AliESDEvent or AliAODEvent !");
180     return ;
181   }
182   fEvent = (AliVEvent*) event;
183
184 }
185
186 //______________________________________________________
187 Bool_t AliHFEextraCuts::IsSelected(TObject *o){
188   //
189   // Steering function for the track selection
190   //
191   TClass *type = o->IsA();
192   AliDebug(2, Form("Object type %s", type->GetName()));
193   if((type == AliESDtrack::Class()) || (type == AliAODTrack::Class())){
194     return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
195   }
196   return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
197 }
198
199 //______________________________________________________
200 Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
201   //
202   // Checks cuts on reconstructed tracks
203   // returns true if track is selected
204   // QA histograms are filled before track selection and for
205   // selected tracks after track selection
206   //
207   AliDebug(1, "Called");
208   ULong64_t survivedCut = 0;    // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
209   if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
210   // Apply cuts
211   Float_t impactR = -999.;
212   Float_t impactZ = -999.;
213   Double_t hfeimpactR, hfeimpactnsigmaR;
214   Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
215   Double_t maximpactRcut; 
216   GetImpactParameters(track, impactR, impactZ);
217   if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
218     // Protection for PbPb
219     GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
220     GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
221   }
222   Int_t  nclsITS = GetITSNbOfcls(track);
223   UInt_t nclsTPC = GetTPCncls(track);
224   UInt_t nclsTPCPID = track->GetTPCsignalN();
225   // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
226   Float_t fractionSharedClustersTPC = GetTPCsharedClustersRatio(track);
227   Double_t ratioTPC = GetTPCclusterRatio(track);
228   UChar_t trdTracklets;
229   trdTracklets = track->GetTRDntrackletsPID();
230   UChar_t itsPixel = track->GetITSClusterMap();
231   Int_t status1 = GetITSstatus(track, 0);
232   Int_t status2 = GetITSstatus(track, 1);
233   Bool_t statusL0 = CheckITSstatus(status1);
234   Bool_t statusL1 = CheckITSstatus(status2);
235   Double_t tofsignalDx = 0.0;
236   Double_t tofsignalDz = 0.0;
237   GetTOFsignalDxDz(track,tofsignalDx,tofsignalDz);
238
239   if(TESTBIT(fRequirements, kTPCfractionShared)) {
240     // cut on max fraction of shared TPC clusters
241     if(TMath::Abs(fractionSharedClustersTPC) < fFractionTPCShared) SETBIT(survivedCut, kTPCfractionShared);    
242   }
243   if(TESTBIT(fRequirements, kMinImpactParamR)){
244     // cut on min. Impact Parameter in Radial direction
245     if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
246   }
247   if(TESTBIT(fRequirements, kMinImpactParamZ)){
248     // cut on min. Impact Parameter in Z direction
249     if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
250   }
251   if(TESTBIT(fRequirements, kMaxImpactParamR)){
252     // cut on max. Impact Parameter in Radial direction
253     if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
254   }
255   if(TESTBIT(fRequirements, kMaxImpactParameterRpar)) {
256     // HFE Impact parameter cut
257     GetMaxImpactParameterCutR(track,maximpactRcut);
258     if(TMath::Abs(impactR) >= maximpactRcut) SETBIT(fRequirements, kMaxImpactParameterRpar);
259   }
260   if(TESTBIT(fRequirements, kMaxImpactParamZ)){
261     // cut on max. Impact Parameter in Z direction
262     if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
263   }
264   if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
265     // cut on min. HFE Impact Parameter in Radial direction
266     if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
267   }
268   if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
269     // cut on max. HFE Impact Parameter n sigma in Radial direction
270     if(fAbsHFEImpactParamNsigmaR) {
271 //      if((TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) && (TMath::Abs(hfeimpactnsigmaR) < 8)) { //mj debug
272       if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) {
273         SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
274 //        printf("0: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
275       }
276     }
277     else {
278       if(hfeimpactnsigmaRcut>0){
279         if(hfeimpactnsigmaR >= hfeimpactnsigmaRcut) {
280           SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
281           //printf("1: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
282         }
283       }
284       else{
285         if(hfeimpactnsigmaR <= hfeimpactnsigmaRcut) {
286           SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
287           //printf("2: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
288         }
289       }
290     }
291   }
292   if(TESTBIT(fRequirements, kClusterRatioTPC)){
293     // cut on min ratio of found TPC clusters vs findable TPC clusters
294     if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
295   }
296   if(TESTBIT(fRequirements, kMinTrackletsTRD)){
297     // cut on minimum number of TRD tracklets
298     AliDebug(1, Form("Min TRD cut: [%d|%d], Require exact number [%s]\n", fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
299     if(fTRDtrackletsExact){
300       if(trdTracklets == fMinTrackletsTRD) {
301         SETBIT(survivedCut, kMinTrackletsTRD);
302         AliDebug(1, "Track Selected");
303       }
304     }else{
305       if(trdTracklets >= fMinTrackletsTRD){ 
306         SETBIT(survivedCut, kMinTrackletsTRD);
307         AliDebug(1, "Track Selected");
308       }
309       //printf("Min number of tracklets %d\n",fMinTrackletsTRD);
310     }
311   }
312   
313   if(TESTBIT(fRequirements, kMinNbITScls)){
314     // cut on minimum number of ITS clusters
315     //printf(Form("Min ITS clusters: [%d|%d]\n", (Int_t)fMinNbITScls, nclsITS));
316     AliDebug(1, Form("Min ITS clusters: [%d|%d]\n", fMinNbITScls, nclsITS));
317     if(nclsITS >= ((Int_t)fMinNbITScls)) SETBIT(survivedCut, kMinNbITScls);
318   }
319   
320   if(TESTBIT(fRequirements, kMinNClustersTPC)){
321     // cut on minimum number of TPC tracklets
322     //printf(Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
323     AliDebug(1, Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
324     if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
325   }
326   if(TESTBIT(fRequirements, kMinNClustersTPCPID)){
327     AliDebug(1, Form("Min TPC PID cut: [%d|%d]\n", fMinNClustersTPCPID, nclsTPCPID));
328     if(nclsTPCPID >= fMinNClustersTPCPID) SETBIT(survivedCut, kMinNClustersTPCPID);
329   }
330   if(TESTBIT(fRequirements, kDriftITS)){
331     //printf("Require drift\n");
332     switch(fDriftITS){
333     case kFirstD:
334       if(TESTBIT(itsPixel, 2)) SETBIT(survivedCut, kDriftITS);
335       break;
336     default: 
337       SETBIT(survivedCut, kDriftITS);
338       break;
339   }
340   }
341   if(TESTBIT(fRequirements, kPixelITS)){
342     // cut on ITS pixel layers
343     AliDebug(1, "ITS cluster Map: ");
344     //PrintBitMap(itsPixel);
345     switch(fPixelITS){
346       case kFirst:
347         AliDebug(2, "First");
348         //printf("First\n");
349               if(TESTBIT(itsPixel, 0)) 
350                 SETBIT(survivedCut, kPixelITS);
351         else
352                 if(fCheck && !statusL0)
353                   SETBIT(survivedCut, kPixelITS);
354                     break;
355       case kSecond: 
356         //printf("Second\n");
357         AliDebug(2, "Second");
358               if(TESTBIT(itsPixel, 1))
359                 SETBIT(survivedCut, kPixelITS);
360         else
361                 if(fCheck && !statusL1)
362                   SETBIT(survivedCut, kPixelITS);
363                     break;
364       case kBoth: 
365         //printf("Both\n");
366         AliDebug(2, "Both");
367               if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
368                       SETBIT(survivedCut, kPixelITS);
369               else
370           if(fCheck && !(statusL0 && statusL1)) 
371                         SETBIT(survivedCut, kPixelITS);
372               break;
373       case kAny: 
374         //printf("Any\n");
375         AliDebug(2, "Any");
376               if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
377                 SETBIT(survivedCut, kPixelITS);
378         else
379                 if(fCheck && !(statusL0 || statusL1))
380                     SETBIT(survivedCut, kPixelITS);
381                     break;
382       case kExclusiveSecond:
383         AliDebug(2, "Exlusive second");
384         if(fCheck){ // Cut out tracks which pass a dead ITS Layer 0
385           if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0) && statusL0)
386             SETBIT(survivedCut, kPixelITS);
387         } else {
388           if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0))
389             SETBIT(survivedCut, kPixelITS);
390         }
391         break;
392       case kNone:
393         // No cut applied, set as survived
394         SETBIT(survivedCut, kPixelITS);
395         break;
396       default: 
397         // default, defined as no cut applied 
398         AliDebug(2, "None");
399         SETBIT(survivedCut, kPixelITS);
400         break;
401     }
402     AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
403   }
404
405   if(TESTBIT(fRequirements, kTOFPID)){
406     // cut on TOF pid
407     if(track->GetStatus() & AliESDtrack::kTOFpid) SETBIT(survivedCut, kTOFPID);
408   }
409
410   if(TESTBIT(fRequirements, kTOFmismatch)){
411     // cut on TOF mismatch
412     if(!(track->GetStatus() & AliESDtrack::kTOFmismatch)) SETBIT(survivedCut, kTOFmismatch);
413   }
414
415   if(TESTBIT(fRequirements, kTPCPIDCleanUp)){
416       // cut on TPC PID cleanup
417       Bool_t fBitsAboveThreshold=GetTPCCountSharedMapBitsAboveThreshold(track);
418       if((fBitsAboveThreshold==0)&&(nclsTPCPID>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
419   }
420
421   if(TESTBIT(fRequirements, kEMCALmatch)){
422     if(track->GetStatus() && AliESDtrack::kEMCALmatch) SETBIT(survivedCut, kEMCALmatch);
423   }
424
425   if(TESTBIT(fRequirements, kRejectKinkDaughter)){
426     //printf("test daughter\n");
427     if(!IsKinkDaughter(track)) SETBIT(survivedCut, kRejectKinkDaughter);
428     //else printf("Found kink daughter\n");
429   }
430
431   if(TESTBIT(fRequirements, kRejectKinkMother)){
432     //printf("test mother\n");
433     if(!IsKinkMother(track)) SETBIT(survivedCut, kRejectKinkMother);
434     //else printf("Found kink mother\n");
435   }
436   if(TESTBIT(fRequirements, kTOFsignalDxy)){
437     // cut on TOF matching cluster
438     if((TMath::Abs(tofsignalDx) <= fTOFsignalDx) && (TMath::Abs(tofsignalDz) <= fTOFsignalDz)) SETBIT(survivedCut, kTOFsignalDxy);
439   }
440   
441   if(fRequirements == survivedCut){
442     //
443     // Track selected
444     //
445     AliDebug(2, "Track Survived cuts\n");
446     if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
447     return kTRUE;
448   }
449   AliDebug(2, "Track cut");
450   if(IsQAOn()) FillCutCorrelation(survivedCut);
451   return kFALSE;
452 }
453
454 //______________________________________________________
455 Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
456   //
457   // Checks cuts on Monte Carlo tracks
458   // returns true if track is selected
459   // QA histograms are filled before track selection and for
460   // selected tracks after track selection
461   //
462   return kTRUE; // not yet implemented
463 }
464
465 //______________________________________________________
466 void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
467   //
468   // Fill the QA histograms for ESD tracks
469   // Function can be called before cuts or after cut application (second argument)
470   //
471   Float_t impactR, impactZ;
472   GetImpactParameters(track, impactR, impactZ);
473   TH1 *htmp = NULL;
474   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * fgkNQAhistos)))) htmp->Fill(impactR);
475   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * fgkNQAhistos)))) htmp->Fill(impactZ);
476   // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
477   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * fgkNQAhistos)))) htmp->Fill(GetTPCclusterRatio(track));
478   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * fgkNQAhistos)))) htmp->Fill(track->GetTRDntrackletsPID());
479   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * fgkNQAhistos)))) htmp->Fill(GetTPCncls(track));
480   UChar_t itsPixel = track->GetITSClusterMap();
481   TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * fgkNQAhistos));
482   //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
483   if(pixelHist){
484     Double_t firstEntry = 0.5;
485     if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
486       pixelHist->Fill(firstEntry + 3);
487     else{
488       if(itsPixel & BIT(0)){
489         pixelHist->Fill(firstEntry);
490         if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
491         else pixelHist->Fill(firstEntry + 4);
492       }
493       if(itsPixel & BIT(1)){
494         pixelHist->Fill(firstEntry + 1);
495         if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
496       }
497     }
498   }
499   // Fill histogram with the status bits
500   TH1 *hStatusBits = dynamic_cast<TH1 *>(fQAlist->At(6 + when * fgkNQAhistos));
501   if(hStatusBits) {
502     hStatusBits->Fill(0);  // Fill first bin with all tracks
503     if(track->GetStatus() && AliESDtrack::kTOFpid) hStatusBits->Fill(1);
504     if(!(track->GetStatus() && AliESDtrack::kTOFmismatch)) hStatusBits->Fill(2);
505     if(track->GetStatus() && AliESDtrack::kEMCALmatch) hStatusBits->Fill(3);
506     if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
507   }
508   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(track->GetTPCsignalN());
509 }
510
511 // //______________________________________________________
512 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
513 //   //
514 //   // Fill the QA histograms for MC tracks
515 //   // Function can be called before cuts or after cut application (second argument)
516 //   // Not yet implemented
517 //   //
518 // }
519
520 //______________________________________________________
521 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
522   //
523   // Fill cut correlation histograms for tracks that didn't pass cuts
524   //
525   TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * fgkNQAhistos));
526   if(correlation){
527     for(Int_t icut = 0; icut < kNcuts; icut++){
528       if(!TESTBIT(fRequirements, icut)) continue;
529       for(Int_t jcut = icut; jcut < kNcuts; jcut++){
530         if(!TESTBIT(fRequirements, jcut)) continue;
531         if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
532                 correlation->Fill(icut, jcut);
533       }
534     }
535   }
536 }
537
538 //______________________________________________________
539 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
540   //
541   // Add QA histograms
542   // For each cut a histogram before and after track cut is created
543   // Histos before respectively after cut are stored in different lists
544   // Additionally a histogram with the cut correlation is created and stored
545   // in the top directory
546   //
547
548   TH1 *histo1D = 0x0;
549   TH2 *histo2D = 0x0;
550   TString cutstr[2] = {"before", "after"};
551
552   if(!fQAlist) fQAlist = new TList;  // for internal representation, not owner
553   for(Int_t icond = 0; icond < 2; icond++){
554     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * fgkNQAhistos);
555     fQAlist->AddAt(histo1D, 0 + icond * fgkNQAhistos);
556     histo1D->GetXaxis()->SetTitle("Impact Parameter");
557     histo1D->GetYaxis()->SetTitle("Number of Tracks");
558     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * fgkNQAhistos);
559     fQAlist->AddAt(histo1D, 1 + icond * fgkNQAhistos);
560     histo1D->GetXaxis()->SetTitle("Impact Parameter");
561     histo1D->GetYaxis()->SetTitle("Number of Tracks");
562     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1)), 2 + icond * fgkNQAhistos);
563     fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
564     histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
565     histo1D->GetYaxis()->SetTitle("Number of Tracks");
566     qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * fgkNQAhistos);
567     fQAlist->AddAt(histo1D, 3 + icond * fgkNQAhistos);
568     histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
569     histo1D->GetYaxis()->SetTitle("Number of Tracks");
570     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * fgkNQAhistos);
571     fQAlist->AddAt(histo1D, 4 + icond * fgkNQAhistos);
572     histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
573     histo1D->GetYaxis()->SetTitle("Number of Tracks");
574     qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * fgkNQAhistos);
575     fQAlist->AddAt(histo1D, 5 + icond * fgkNQAhistos);
576     histo1D->GetXaxis()->SetTitle("ITS Pixel");
577     histo1D->GetYaxis()->SetTitle("Number of Tracks");
578     Int_t first = histo1D->GetXaxis()->GetFirst();
579     TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
580     for(Int_t ilabel = 0; ilabel < 6; ilabel++)
581       histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
582     qaList->AddAt((histo1D = new TH1F(Form("%s_trackStatus%s",GetName(),cutstr[icond].Data()), "Track Status", 5, 0, 5)), 6 + icond * fgkNQAhistos);
583     fQAlist->AddAt(histo1D, 6 + icond * fgkNQAhistos);
584     histo1D->GetXaxis()->SetTitle("Track Status Bit");
585     histo1D->GetYaxis()->SetTitle("Number of tracks");
586     TString tsnames[5] = {"All", "TOFPID", "No TOFmismatch", "EMCALmatch","No TPC shared bit"};
587     for(Int_t istat = 0; istat < 5; istat++) histo1D->GetXaxis()->SetBinLabel(istat + 1, tsnames[istat].Data());
588     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcdEdxClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters for dEdx calculation", 161, 0, 160)), 7 + icond * fgkNQAhistos);
589     fQAlist->AddAt(histo1D, 7 + icond * fgkNQAhistos);
590     histo1D->GetXaxis()->SetTitle("Number of TPC clusters for dEdx calculation");
591     histo1D->GetYaxis()->SetTitle("counts");
592   }
593   // Add cut correlation
594   qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * fgkNQAhistos);
595   fQAlist->AddAt(histo2D, 2 * fgkNQAhistos);
596   TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel", "kMinHFEImpactParamR", "kMinHFEImpactParamNsigmaR", "TPC Number of clusters" "Fraction Shared TPC clusters", "TOFPID", "No TOFmismatch", "EMCALmatch", "ImpactParam"};
597   Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
598   for(Int_t icut = 0; icut < kNcuts; icut++){
599     histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
600     histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
601   }
602 }
603
604 //______________________________________________________
605 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
606   for(Int_t ibit = 32; ibit--; )
607     printf("%d", bitmap & BIT(ibit) ? 1 : 0);
608   printf("\n");
609 }
610
611 //______________________________________________________
612 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
613   //
614   // Check whether ITS area is dead
615   //
616   Bool_t status;
617   switch(itsStatus){
618     case 2: status = kFALSE; break;
619     case 3: status = kFALSE; break;
620     case 7: status = kFALSE; break;
621     default: status = kTRUE;
622   }
623   return status;
624 }
625
626 //______________________________________________________
627 Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
628         //
629         // Check ITS layer status
630         //
631         Int_t status = 0;
632         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
633                 Int_t det;
634                 Float_t xloc, zloc;
635                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
636                 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
637         }
638         return status;
639 }
640
641
642 //______________________________________________________
643 Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track){
644   //
645   // Checks if number of shared bits is above 1
646   //
647   Int_t nsharebit = 1;
648   const TBits *shared;
649   if(track->IsA() == AliESDtrack::Class())
650     shared = &((static_cast<AliESDtrack *>(track))->GetTPCSharedMap());
651   else if(track->IsA() == AliAODTrack::Class())
652     shared = &((static_cast<AliAODTrack *>(track))->GetTPCSharedMap());
653    else 
654     return 0;
655
656         if(shared->CountBits() >= 2) nsharebit=1;
657   else nsharebit=0;
658  
659   return nsharebit;
660 }
661
662 //______________________________________________________
663 UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
664         //
665         // Get Number of findable clusters in the TPC
666         //
667         Int_t nClusters = 0; // in case no Information available consider all clusters findable
668   TClass *type = track->IsA();
669   if(TESTBIT(fTPCclusterDef, kFoundIter1)){
670                 if(type == AliESDtrack::Class()){
671       AliDebug(2, ("Using def kFoundIter1"));
672       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
673       nClusters = esdtrack->GetTPCNclsIter1();
674     } else {
675       AliDebug(2, ("Number of clusters in iteration 1 not available for AOD tracks"));
676     }
677   } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
678     AliDebug(2, ("Using def kCrossedRows"));
679     nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
680   } else{
681     AliDebug(2, ("Using def kFound"));
682           nClusters = track->GetTPCNcls();
683   }
684         return nClusters;
685 }
686
687 //______________________________________________________
688 Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
689   //
690   // Get fraction of shared TPC clusters
691   //
692   Float_t fracClustersTPCShared = 0.0; 
693   Int_t nClustersTPC = track->GetTPCNcls();
694   Int_t nClustersTPCShared = 0;
695   TClass *type = track->IsA();
696   if(type == AliESDtrack::Class()){
697     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
698     nClustersTPCShared = esdtrack->GetTPCnclsS();
699   } else if(type == AliESDtrack::Class()){
700     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
701     const TBits &shared = aodtrack->GetTPCSharedMap();
702     nClustersTPCShared = shared.CountBits(0) - shared.CountBits(159);
703   }
704   if (nClustersTPC!=0) {
705            fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
706   }
707   return fracClustersTPCShared;
708 }
709
710 //______________________________________________________
711 Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
712   //
713   // Get Ratio of found / findable clusters for different definitions
714   // Only implemented for ESD tracks
715   //
716   Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
717   TClass *type = track->IsA();
718   if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
719     AliDebug(2, "Using ratio def kCROverFindable");
720     clusterRatio = track->GetTPCNclsF() ? track->GetTPCClusterInfo(2,1)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // crossed rows/findable
721   } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
722     AliDebug(2, "Using ratio def kFoundOverCR");
723     clusterRatio = track->GetTPCClusterInfo(2,0); // found/crossed rows
724   } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
725     if(type == AliESDtrack::Class()){
726       AliDebug(2, "Using ratio def kFoundOverFindableIter1");
727       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
728       clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
729     } else {
730       AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
731       clusterRatio = 1.;
732     }
733   } else {
734     AliDebug(2, "Using ratio def kFoundOverFindable");
735     clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // found/findable
736   }
737   return clusterRatio;
738 }
739
740 //______________________________________________________
741 void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
742   //
743   // Get impact parameter
744   //
745   TClass *type = track->IsA();
746   if(type == AliESDtrack::Class()){
747     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
748     esdtrack->GetImpactParameters(radial, z);
749   }
750   else if(type == AliAODTrack::Class()){
751     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
752     if(aodtrack){
753       Double_t xyz[3];
754       aodtrack->XYZAtDCA(xyz);
755       //printf("xyz[0] %f, xyz[1] %f, xyz[2] %f\n",xyz[0], xyz[1],xyz[2]);
756       if(xyz[0]==-999. || xyz[1]==-999. ||  xyz[2]==-999.){
757         if(!fEvent) {
758           z = xyz[2];
759           radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
760           //printf("No event\n");
761         }
762         else {
763           // Propagate the global track to the DCA.
764           const AliVVertex *vAOD = fEvent->GetPrimaryVertex();
765           Double_t PosAtDCA[2] = {-999,-999};
766           Double_t covar[3] = {-999,-999,-999};
767           AliAODTrack *copyaodtrack = new AliAODTrack(*aodtrack);
768           if(copyaodtrack && copyaodtrack->PropagateToDCA(vAOD,fEvent->GetMagneticField(),100.,PosAtDCA,covar)) {
769             z = PosAtDCA[1];
770             radial = PosAtDCA[0];
771             //printf("Propagate z %f and radial %f\n",z,radial);
772           }
773           delete copyaodtrack;
774         }
775       }
776       else {
777         z = xyz[2];
778         radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
779       }
780     }
781   }
782 }
783 //______________________________________________________
784 Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
785   //
786   // Is kink Daughter
787   //
788   TClass *type = track->IsA();
789   if(type == AliESDtrack::Class()){
790     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
791     if(!esdtrack) return kFALSE;
792     if(esdtrack->GetKinkIndex(0)<=0) return kFALSE;
793     else return kTRUE;
794
795   }
796   else if(type == AliAODTrack::Class()){
797     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
798     if(aodtrack){
799       AliAODVertex *aodvertex = aodtrack->GetProdVertex();
800       if(!aodvertex) return kFALSE;
801       if(aodvertex->GetType()==AliAODVertex::kKink) return kTRUE;
802       else return kFALSE;
803     }
804     else return kFALSE;
805   }
806
807   return kFALSE;
808 }
809 //______________________________________________________
810 Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
811   //
812   // Is kink Mother: only for ESD since need to loop over vertices for AOD
813   //
814   //
815
816   TClass *type = track->IsA();
817   if(type == AliESDtrack::Class()){
818     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
819     if(!esdtrack) return kFALSE;
820     if(esdtrack->GetKinkIndex(0)!=0) return kTRUE;
821     else return kFALSE;
822   }
823
824   return kFALSE;
825
826 }
827 //______________________________________________________
828 Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
829   //
830   // Get ITS nb of clusters
831   //
832   TClass *type = track->IsA();
833   if(type == AliESDtrack::Class()){
834     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
835     return esdtrack->GetITSclusters(0);
836
837   }
838   else if(type == AliAODTrack::Class()){
839     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
840     if(aodtrack){
841       return  aodtrack->GetITSNcls();
842     }
843   }
844
845   return 0;
846 }
847 //______________________________________________________
848 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
849         //
850         // Get HFE impact parameter (with recalculated primary vertex)
851         //
852         dcaxy=0;
853         dcansigmaxy=0;
854   if(!fEvent){
855     AliDebug(1, "No Input event available\n");
856     return;
857   }
858   const Double_t kBeampiperadius=3.;
859   TString type = track->IsA()->GetName();
860   Double_t dcaD[2]={-999.,-999.};
861   Double_t covD[3]={-999.,-999.,-999.};
862   Float_t dcaF[2]={-999.,-999.};
863   Float_t covF[3]={-999.,-999.,-999.};
864   
865   AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
866   if(!esdevent) {
867     AliDebug(1, "No esd event available\n");
868     return;
869   }
870   const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
871   if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
872    // recalculate primary vertex for peri. and pp
873    AliVertexerTracks vertexer(fEvent->GetMagneticField());
874    vertexer.SetITSMode();
875    vertexer.SetMinClusters(4);
876          Int_t skipped[2];
877    skipped[0] = track->GetID();
878    vertexer.SetSkipTracks(1,skipped);
879    
880  //----diamond constraint-----------------------------
881    vertexer.SetConstraintOn();
882    Float_t diamondcovxy[3];
883    esdevent->GetDiamondCovXY(diamondcovxy);
884    Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
885    Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
886    AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
887    vertexer.SetVtxStart(diamond);
888    delete diamond; diamond=NULL;
889  //----------------------------------------------------
890     
891    vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
892
893    // Getting the DCA
894    // Propagation always done on a working copy to not disturb the track params of the original track
895    AliESDtrack *esdtrack = NULL;
896    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
897      // Case ESD track: take copy constructor
898      AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
899      if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
900    } else {
901     // Case AOD track: take different constructor
902     esdtrack = new AliESDtrack(track);
903    }
904    if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD)){
905     // protection
906     dcaxy = dcaD[0];
907     if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
908     if(!covD[0]) dcansigmaxy = -49.;
909    }
910    delete esdtrack;
911    delete vtxESDSkip;
912   } 
913   else{
914    AliESDtrack *esdtrack = NULL;
915    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
916     // Case ESD track: take copy constructor
917     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
918     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
919    } else {
920     // Case AOD track: take different constructor
921     esdtrack = new AliESDtrack(track);
922    }
923    if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF); 
924    dcaxy = dcaF[0];
925    if(covF[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covF[0]);
926    if(!covF[0]) dcansigmaxy = -49.;
927    delete esdtrack;
928   }
929 }
930
931 //______________________________________________________
932 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t dcaD[2], Double_t covD[3]){
933         //
934         // Get HFE impact parameter (with recalculated primary vertex)
935         //
936   if(!fEvent){
937     AliDebug(1, "No Input event available\n");
938     return;
939   }
940   const Double_t kBeampiperadius=3.;
941   TString type = track->IsA()->GetName();
942   Float_t dcaF[2]={-999.,-999.};
943   Float_t covF[3]={-999.,-999.,-999.};
944   
945   AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
946   if(!esdevent) {
947     AliDebug(1, "No esd event available\n");
948     return;
949   }
950   const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
951   if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
952    // recalculate primary vertex for peri. and pp
953    AliVertexerTracks vertexer(fEvent->GetMagneticField());
954    vertexer.SetITSMode();
955    vertexer.SetMinClusters(4);
956          Int_t skipped[2];
957    skipped[0] = track->GetID();
958    vertexer.SetSkipTracks(1,skipped);
959    
960  //----diamond constraint-----------------------------
961    vertexer.SetConstraintOn();
962    Float_t diamondcovxy[3];
963    esdevent->GetDiamondCovXY(diamondcovxy);
964    Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
965    Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
966    AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
967    vertexer.SetVtxStart(diamond);
968    delete diamond; diamond=NULL;
969  //----------------------------------------------------
970     
971    vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
972
973    // Getting the DCA
974    // Propagation always done on a working copy to not disturb the track params of the original track
975    AliESDtrack *esdtrack = NULL;
976    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
977      // Case ESD track: take copy constructor
978      AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
979      if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
980    } else {
981     // Case AOD track: take different constructor
982     esdtrack = new AliESDtrack(track);
983    }
984    if(esdtrack) esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD);
985    delete esdtrack;
986    delete vtxESDSkip;
987   } 
988   else{
989    AliESDtrack *esdtrack = NULL;
990    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
991     // Case ESD track: take copy constructor
992     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
993     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
994    } else {
995     // Case AOD track: take different constructor
996     esdtrack = new AliESDtrack(track);
997    }
998    if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF); 
999    dcaD[0] = dcaF[0];
1000    dcaD[1] = dcaF[1];
1001    covD[0] = covF[0];
1002    covD[1] = covF[1];
1003    covD[2] = covF[2];
1004    delete esdtrack;
1005   }
1006 }
1007
1008 //______________________________________________________
1009 void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
1010         //
1011         // Get HFE impact parameter cut(pt dependent)
1012         //
1013   
1014   TString type = track->IsA()->GetName();
1015   if(!type.CompareTo("AliESDtrack")){
1016     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
1017     if(!esdtrack) return;
1018     Double_t pt = esdtrack->Pt();       
1019     hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt);  // abs R cut
1020     hfeimpactnsigmaRcut = fIPcutParam[3];                                  // sigma cut
1021   }
1022 }
1023 //______________________________________________________
1024 void AliHFEextraCuts::GetMaxImpactParameterCutR(AliVTrack *track, Double_t &maximpactRcut){
1025         //
1026         // Get max impact parameter cut r (pt dependent)
1027         //
1028   
1029   TString type = track->IsA()->GetName();
1030   if(!type.CompareTo("AliESDtrack")){
1031     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
1032     if(!esdtrack) return;
1033     Double_t pt = esdtrack->Pt();       
1034     if(pt > 0.15) {
1035       maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01);  // abs R cut
1036     }
1037     else maximpactRcut = 9999999999.0;
1038   }
1039 }
1040 //______________________________________________________
1041 void AliHFEextraCuts::GetTOFsignalDxDz(AliVTrack *track, Double_t &tofsignalDx, Double_t &tofsignalDz){
1042   //
1043   // TOF matching 
1044   //
1045   
1046   TString type = track->IsA()->GetName();
1047   if(!type.CompareTo("AliESDtrack")){
1048     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
1049     if(!esdtrack) return;
1050     tofsignalDx = esdtrack->GetTOFsignalDx();
1051     tofsignalDz = esdtrack->GetTOFsignalDz();
1052   }
1053
1054 }