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