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