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, kTPCPIDCleanUp)){
 467       // cut on TPC PID cleanup
 468       Bool_t fBitsAboveThreshold=GetTPCCountSharedMapBitsAboveThreshold(track);
 469       if((fBitsAboveThreshold==0)&&(nclsTPCPID>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
 470   }
 471
 472   if(TESTBIT(fRequirements, kEMCALmatch)){
 473     if(track->GetStatus() && AliESDtrack::kEMCALmatch) SETBIT(survivedCut, kEMCALmatch);
 474   }
 475
 476   if(TESTBIT(fRequirements, kRejectKinkDaughter)){
 477     //printf("test daughter\n");
 478     if(!IsKinkDaughter(track)) SETBIT(survivedCut, kRejectKinkDaughter);
 479     //else printf("Found kink daughter\n");
 480   }
 481
 482   if(TESTBIT(fRequirements, kRejectKinkMother)){
 483     //printf("test mother\n");
 484     if(!IsKinkMother(track)) SETBIT(survivedCut, kRejectKinkMother);
 485     //else printf("Found kink mother\n");
 486   }
 487   if(TESTBIT(fRequirements, kTOFsignalDxy)){
 488     // cut on TOF matching cluster
 489     if((TMath::Abs(tofsignalDx) <= fTOFsignalDx) && (TMath::Abs(tofsignalDz) <= fTOFsignalDz)) SETBIT(survivedCut, kTOFsignalDxy);
 490   }
 491   if(TESTBIT(fRequirements, kITSpattern)){
 492     // cut on ITS pattern (every layer with a working ITS module must have an ITS cluster)
 493     if(CheckITSpattern(track)) SETBIT(survivedCut, kITSpattern);
 494   }
 495   if(TESTBIT(fRequirements, kAODFilterBit)){
 496     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 497     if(aodtrack){
 498       Int_t aodfilter = 1 << fAODFilterBit;
 499       if(aodtrack->TestFilterBit(aodfilter)) SETBIT(survivedCut, kAODFilterBit);
 500     }
 501   }
 502
 503   if(fRequirements == survivedCut){
 504     //
 505     // Track selected
 506     //
 507     AliDebug(2, Form("%s: Track Survived cuts\n", GetName()));
 508     if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
 509     return kTRUE;
 510   }
 511   AliDebug(2, Form("%s: Track cut", GetName()));
 512   if(IsQAOn()) FillCutCorrelation(survivedCut);
 513   return kFALSE;
 514 }
 515
 516 //______________________________________________________
 517 Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
 518   //
 519   // Checks cuts on Monte Carlo tracks
 520   // returns true if track is selected
 521   // QA histograms are filled before track selection and for
 522   // selected tracks after track selection
 523   //
 524   return kTRUE; // not yet implemented
 525 }
 526
 527 //______________________________________________________
 528 void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
 529   //
 530   // Fill the QA histograms for ESD tracks
 531   // Function can be called before cuts or after cut application (second argument)
 532   //
 533   Float_t impactR, impactZ;
 534   GetImpactParameters(track, impactR, impactZ);
 535   TH1 *htmp = NULL;
 536   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * fgkNQAhistos)))) htmp->Fill(impactR);
 537   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * fgkNQAhistos)))) htmp->Fill(impactZ);
 538   // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 539   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * fgkNQAhistos)))) htmp->Fill(GetTPCclusterRatio(track));
 540   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * fgkNQAhistos)))) htmp->Fill(track->GetTRDntrackletsPID());
 541   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * fgkNQAhistos)))) htmp->Fill(GetTPCncls(track));
 542   UChar_t itsPixel = track->GetITSClusterMap();
 543   TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * fgkNQAhistos));
 544   //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
 545   if(pixelHist){
 546     Double_t firstEntry = 0.5;
 547     if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
 548       pixelHist->Fill(firstEntry + 3);
 549     else{
 550       if(itsPixel & BIT(0)){
 551         pixelHist->Fill(firstEntry);
 552         if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
 553         else pixelHist->Fill(firstEntry + 4);
 554       }
 555       if(itsPixel & BIT(1)){
 556         pixelHist->Fill(firstEntry + 1);
 557         if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
 558       }
 559     }
 560   }
 561   // Fill histogram with the status bits
 562   TH1 *hStatusBits = dynamic_cast<TH1 *>(fQAlist->At(6 + when * fgkNQAhistos));
 563   if(hStatusBits) {
 564     hStatusBits->Fill(0);  // Fill first bin with all tracks
 565     if(track->GetStatus() && AliESDtrack::kTOFpid) hStatusBits->Fill(1);
 566     if(!(track->GetStatus() && AliESDtrack::kTOFmismatch)) hStatusBits->Fill(2);
 567     if(track->GetStatus() && AliESDtrack::kEMCALmatch) hStatusBits->Fill(3);
 568     if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
 569   }
 570   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(track->GetTPCsignalN());
 571
 572   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(8 + when * fgkNQAhistos)))) htmp->Fill(GetTRDchi(track));
 573
 574   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(9 + when * fgkNQAhistos)))) htmp->Fill(GetITSNbOfcls(track));
 575 }
 576
 577 // //______________________________________________________
 578 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
 579 //   //
 580 //   // Fill the QA histograms for MC tracks
 581 //   // Function can be called before cuts or after cut application (second argument)
 582 //   // Not yet implemented
 583 //   //
 584 // }
 585
 586 //______________________________________________________
 587 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
 588   //
 589   // Fill cut correlation histograms for tracks that didn't pass cuts
 590   //
 591   TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * fgkNQAhistos));
 592   if(correlation){
 593     for(Int_t icut = 0; icut < kNcuts; icut++){
 594       if(!TESTBIT(fRequirements, icut)) continue;
 595       for(Int_t jcut = icut; jcut < kNcuts; jcut++){
 596         if(!TESTBIT(fRequirements, jcut)) continue;
 597         if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
 598                 correlation->Fill(icut, jcut);
 599       }
 600     }
 601   }
 602 }
 603
 604 //______________________________________________________
 605 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
 606   //
 607   // Add QA histograms
 608   // For each cut a histogram before and after track cut is created
 609   // Histos before respectively after cut are stored in different lists
 610   // Additionally a histogram with the cut correlation is created and stored
 611   // in the top directory
 612   //
 613
 614   TH1 *histo1D = 0x0;
 615   TH2 *histo2D = 0x0;
 616   TString cutstr[2] = {"before", "after"};
 617
 618   if(!fQAlist) fQAlist = new TList;  // for internal representation, not owner
 619   for(Int_t icond = 0; icond < 2; icond++){
 620     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * fgkNQAhistos);
 621     fQAlist->AddAt(histo1D, 0 + icond * fgkNQAhistos);
 622     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 623     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 624     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * fgkNQAhistos);
 625     fQAlist->AddAt(histo1D, 1 + icond * fgkNQAhistos);
 626     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 627     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 628     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1.3)), 2 + icond * fgkNQAhistos);
 629     fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
 630     histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
 631     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 632     qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * fgkNQAhistos);
 633     fQAlist->AddAt(histo1D, 3 + icond * fgkNQAhistos);
 634     histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
 635     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 636     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * fgkNQAhistos);
 637     fQAlist->AddAt(histo1D, 4 + icond * fgkNQAhistos);
 638     histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
 639     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 640     qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * fgkNQAhistos);
 641     fQAlist->AddAt(histo1D, 5 + icond * fgkNQAhistos);
 642     histo1D->GetXaxis()->SetTitle("ITS Pixel");
 643     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 644     Int_t first = histo1D->GetXaxis()->GetFirst();
 645     TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
 646     for(Int_t ilabel = 0; ilabel < 6; ilabel++)
 647       histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
 648     qaList->AddAt((histo1D = new TH1F(Form("%s_trackStatus%s",GetName(),cutstr[icond].Data()), "Track Status", 5, 0, 5)), 6 + icond * fgkNQAhistos);
 649     fQAlist->AddAt(histo1D, 6 + icond * fgkNQAhistos);
 650     histo1D->GetXaxis()->SetTitle("Track Status Bit");
 651     histo1D->GetYaxis()->SetTitle("Number of tracks");
 652     TString tsnames[5] = {"All", "TOFPID", "No TOFmismatch", "EMCALmatch","No TPC shared bit"};
 653     for(Int_t istat = 0; istat < 5; istat++) histo1D->GetXaxis()->SetBinLabel(istat + 1, tsnames[istat].Data());
 654     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);
 655     fQAlist->AddAt(histo1D, 7 + icond * fgkNQAhistos);
 656     histo1D->GetXaxis()->SetTitle("Number of TPC clusters for dEdx calculation");
 657     histo1D->GetYaxis()->SetTitle("counts");
 658     qaList->AddAt((histo1D = new TH1F(Form("%s_trdchi2perTracklet%s",GetName(),cutstr[icond].Data()), "chi2 per TRD tracklet", 100, 0, 10)), 8 + icond * fgkNQAhistos);
 659     fQAlist->AddAt(histo1D, 8 + icond * fgkNQAhistos);
 660     histo1D->GetXaxis()->SetTitle("Chi2 per TRD Tracklet");
 661     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 662     qaList->AddAt((histo1D = new TH1F(Form("%s_ITScls%s",GetName(),cutstr[icond].Data()), "ITScls", 6, 0, 6)), 9 + icond * fgkNQAhistos);
 663     fQAlist->AddAt(histo1D, 9 + icond * fgkNQAhistos);
 664     histo1D->GetXaxis()->SetTitle("ITScls");
 665     histo1D->GetYaxis()->SetTitle("Number of ITS cls");
 666
 667   }
 668   // Add cut correlation
 669   qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * fgkNQAhistos);
 670   fQAlist->AddAt(histo2D, 2 * fgkNQAhistos);
 671   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"};
 672   Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
 673   for(Int_t icut = 0; icut < kNcuts; icut++){
 674     histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
 675     histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
 676   }
 677 }
 678
 679 //______________________________________________________
 680 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
 681   for(Int_t ibit = 32; ibit--; )
 682     printf("%d", bitmap & BIT(ibit) ? 1 : 0);
 683   printf("\n");
 684 }
 685
 686 //______________________________________________________
 687 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
 688   //
 689   // Check whether ITS area is dead
 690   //
 691   Bool_t status;
 692   switch(itsStatus){
 693     case 2: status = kFALSE; break;
 694     case 3: status = kFALSE; break;
 695     case 7: status = kFALSE; break;
 696     default: status = kTRUE;
 697   }
 698   return status;
 699 }
 700
 701 //______________________________________________________
 702 Int_t AliHFEextraCuts::GetITSstatus(const AliVTrack * const track, Int_t layer) const {
 703         //
 704         // Check ITS layer status
 705         //
 706         Int_t status = 0;
 707         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 708                 Int_t det;
 709                 Float_t xloc, zloc;
 710                 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 711                 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
 712         }
 713         return status;
 714 }
 715
 716
 717 //______________________________________________________
 718 Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track){
 719   //
 720   // Checks if number of shared bits is above 1
 721   //
 722   Int_t nsharebit = 1;
 723   const TBits *shared;
 724   if(track->IsA() == AliESDtrack::Class())
 725     shared = &((static_cast<AliESDtrack *>(track))->GetTPCSharedMap());
 726   else if(track->IsA() == AliAODTrack::Class())
 727     shared = &((static_cast<AliAODTrack *>(track))->GetTPCSharedMap());
 728    else
 729     return 0;
 730
 731         if(shared->CountBits() >= 2) nsharebit=1;
 732   else nsharebit=0;
 733
 734   return nsharebit;
 735 }
 736
 737 //______________________________________________________
 738 UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
 739   //
 740   // Get Number of findable clusters in the TPC
 741   //
 742   Int_t nClusters = 0; // in case no Information available consider all clusters findable
 743   TClass *type = track->IsA();
 744   if(TESTBIT(fTPCclusterDef, kFoundIter1)){
 745     if(type == AliESDtrack::Class()){
 746       AliDebug(2, ("Using def kFoundIter1"));
 747       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 748       nClusters = esdtrack->GetTPCNclsIter1();
 749     } else {
 750       AliDebug(2, ("Number of clusters in iteration 1 not available for AOD tracks"));
 751     }
 752   } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
 753     AliDebug(2, ("Using def kCrossedRows"));
 754     nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
 755   } else if(TESTBIT(fTPCclusterDef, kFound)){
 756     AliDebug(2, ("Using def kFound"));
 757     nClusters = track->GetTPCNcls();
 758   }
 759   else {
 760     AliDebug(2, ("Using def kFoundAll"));
 761     if(type == AliESDtrack::Class()){
 762       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 763       const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
 764       nClusters = clusterTPC.CountBits();
 765     }
 766     else {
 767       AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
 768       const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
 769       nClusters = clusterTPC.CountBits();
 770     }
 771   }
 772   return nClusters;
 773 }
 774
 775 //______________________________________________________
 776 Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
 777   //
 778   // Get fraction of shared TPC clusters
 779   //
 780   Float_t fracClustersTPCShared = 0.0;
 781   Int_t nClustersTPC = track->GetTPCNcls();
 782   Int_t nClustersTPCShared = 0;
 783   TClass *type = track->IsA();
 784   if(type == AliESDtrack::Class()){
 785     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 786     nClustersTPCShared = esdtrack->GetTPCnclsS();
 787   } else if(type == AliAODTrack::Class()){
 788     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
 789     const TBits &shared = aodtrack->GetTPCSharedMap();
 790     nClustersTPCShared = shared.CountBits(0) - shared.CountBits(159);
 791   }
 792   if (nClustersTPC!=0) {
 793     fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
 794   }
 795   return fracClustersTPCShared;
 796 }
 797
 798 //______________________________________________________
 799 Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
 800   //
 801   // Get Ratio of found / findable clusters for different definitions
 802   // Only implemented for ESD tracks
 803   //
 804   Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
 805   TClass *type = track->IsA();
 806   if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
 807     AliDebug(2, "Using ratio def kCROverFindable");
 808     clusterRatio = track->GetTPCNclsF() ? track->GetTPCClusterInfo(2,1)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // crossed rows/findable
 809   } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
 810     AliDebug(2, "Using ratio def kFoundOverCR");
 811     clusterRatio = track->GetTPCClusterInfo(2,0); // found/crossed rows
 812   } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
 813     if(type == AliESDtrack::Class()){
 814       AliDebug(2, "Using ratio def kFoundOverFindableIter1");
 815       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 816       clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
 817     } else {
 818       AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
 819       clusterRatio = 1.;
 820     }
 821   } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindable)) {
 822     AliDebug(2, "Using ratio def kFoundOverFindable");
 823     clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // found/findable
 824   }
 825   else {
 826     AliDebug(2, "Using ratio def kFoundAllOverFindable");
 827     Int_t allclusters = 0;
 828     if(type == AliESDtrack::Class()){
 829       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 830       const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
 831       allclusters = clusterTPC.CountBits();
 832     }
 833     else {
 834       AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
 835       const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
 836       allclusters = clusterTPC.CountBits();
 837     }
 838     clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(allclusters)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // foundall/findable
 839   }
 840   return clusterRatio;
 841
 842 }
 843 //___________________________________________________
 844 void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
 845   //
 846   // Get impact parameter
 847   //
 848
 849   const Double_t kBeampiperadius=3.;
 850   Double_t dcaD[2]={-999.,-999.},
 851            covD[3]={-999.,-999.,-999.};
 852   TClass *type = track->IsA();
 853   if(type == AliESDtrack::Class()){
 854     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 855     esdtrack->GetImpactParameters(radial, z);
 856   }
 857   else if(type == AliAODTrack::Class()){
 858
 859     //case of AOD tracks
 860     AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
 861     if(!aodevent) {
 862       AliDebug(1, "No aod event available\n");
 863       return;
 864     }
 865
 866     //Case ESD track: take copy constructor
 867     AliAODTrack *aodtrack = NULL;
 868     AliAODTrack *tmptrack = dynamic_cast<AliAODTrack *>(track);
 869     if(tmptrack) aodtrack = new AliAODTrack(*tmptrack);
 870
 871     AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
 872     if(!vtxAODSkip) return;
 873     AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
 874     if(etp.PropagateToDCA(vtxAODSkip, aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
 875       radial = dcaD[0];
 876       z = dcaD[1];
 877     }
 878     //if(vtxAODSkip) delete vtxAODSkip;
 879     if(aodtrack) delete aodtrack;
 880   }
 881 }
 882 //______________________________________________________
 883 Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
 884   //
 885   // Is kink Daughter
 886   //
 887   TClass *type = track->IsA();
 888   if(type == AliESDtrack::Class()){
 889     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 890     if(!esdtrack) return kFALSE;
 891     if(esdtrack->GetKinkIndex(0)>0) return kTRUE;
 892     else return kFALSE;
 893
 894   }
 895   else if(type == AliAODTrack::Class()){
 896     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 897     if(aodtrack){
 898       AliAODVertex *aodvertex = aodtrack->GetProdVertex();
 899       if(!aodvertex) return kFALSE;
 900       if(aodvertex->GetType()==AliAODVertex::kKink) return kTRUE;
 901       else return kFALSE;
 902     }
 903     else return kFALSE;
 904   }
 905
 906   return kFALSE;
 907 }
 908 //______________________________________________________
 909 Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
 910   //
 911   // Is kink Mother
 912   //
 913
 914   TClass *type = track->IsA();
 915   if(type == AliESDtrack::Class()){
 916     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 917     if(!esdtrack) return kFALSE;
 918     if(esdtrack->GetKinkIndex(0)!=0) return kTRUE;
 919     else return kFALSE;
 920   } else if(type == AliAODTrack::Class()){
 921     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
 922     for(int ikink = 0; ikink < fNumberKinkMothers; ikink++){
 923       if(fListKinkMothers[ikink] == aodtrack->GetID()) return kTRUE;
 924     }
 925   }
 926
 927   return kFALSE;
 928
 929 }
 930
 931 //______________________________________________________
 932 Float_t AliHFEextraCuts::GetTRDchi(AliVTrack *track){
 933   //
 934   // Get TRDchi2
 935   //
 936   Int_t ntls(0);
 937   TClass *type = track->IsA();
 938   if(type == AliESDtrack::Class()){
 939       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 940       ntls = esdtrack->GetTRDntracklets();
 941       return ntls ? esdtrack->GetTRDchi2()/ntls : -999;
 942   }
 943   else if(type == AliAODTrack::Class()){
 944     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 945     if(aodtrack){
 946       return  999.;
 947     }
 948   }
 949
 950   return 999.;
 951
 952 }
 953
 954 //______________________________________________________
 955 Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
 956   //
 957   // Get ITS nb of clusters
 958   //
 959   TClass *type = track->IsA();
 960   if(type == AliESDtrack::Class()){
 961     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 962     return esdtrack->GetITSclusters(0);
 963
 964   }
 965   else if(type == AliAODTrack::Class()){
 966     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 967     if(aodtrack){
 968       return  aodtrack->GetITSNcls();
 969     }
 970   }
 971
 972   return 0;
 973 }
 974 //______________________________________________________
 975 void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t &dcaxy, Double_t &dcansigmaxy){
 976   //
 977   // Get HFE impact parameter (with recalculated primary vertex)
 978   //
 979   dcaxy=0;
 980   dcansigmaxy=0;
 981   if(!fEvent){
 982     AliDebug(1, "No Input event available\n");
 983     return;
 984   }
 985   TString type = track->IsA()->GetName();
 986   const Double_t kBeampiperadius=3.;
 987   Double_t dcaD[2]={-999.,-999.},
 988            covD[3]={-999.,-999.,-999.};
 989   Bool_t isRecalcVertex(kFALSE);
 990
 991   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 992     //case of ESD tracks
 993     AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
 994     if(!esdevent) {
 995       AliDebug(1, "No esd event available\n");
 996       return;
 997     }
 998
 999     const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
1000     if(!vtxESDSkip) return;
1001
1002     //case ESD track: take copy constructor
1003     const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
1004     if(tmptrack){
1005
1006       if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1007         vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
1008         isRecalcVertex = kTRUE;
1009       }
1010
1011       if(vtxESDSkip){
1012         AliESDtrack esdtrack(*tmptrack);
1013         fMagField = fEvent->GetMagneticField();
1014         if(esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD)){
1015           dcaxy = dcaD[0];
1016           if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
1017         }
1018         if(isRecalcVertex) delete vtxESDSkip;
1019       }
1020     }
1021   }
1022   else {
1023     //case of AOD tracks
1024     AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
1025     if(!aodevent) {
1026       AliDebug(1, "No aod event available\n");
1027       return;
1028     }
1029
1030     AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
1031     if(!vtxAODSkip) return;
1032
1033     //Case ESD track: take copy constructor
1034     const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
1035     if(tmptrack){
1036
1037       if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1038         vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
1039         isRecalcVertex = kTRUE;
1040       }
1041       if(vtxAODSkip){
1042         AliAODTrack aodtrack(*tmptrack);
1043         AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
1044         fMagField = aodevent->GetMagneticField();
1045         if(etp.PropagateToDCA(vtxAODSkip,fMagField, kBeampiperadius, dcaD, covD)) {
1046           dcaxy = dcaD[0];
1047           if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
1048         }
1049         if(isRecalcVertex) delete vtxAODSkip;
1050       }
1051     }
1052   }
1053 }
1054
1055 //______________________________________________________
1056 void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t dcaD[2], Double_t covD[3]){
1057         //
1058         // Get HFE impact parameter (with recalculated primary vertex)
1059         //
1060   if(!fEvent){
1061     AliDebug(1, "No Input event available\n");
1062     return;
1063   }
1064   const Double_t kBeampiperadius=3.;
1065   TString type = track->IsA()->GetName();
1066   Bool_t isRecalcVertex(kFALSE);
1067
1068   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
1069     //case of ESD tracks
1070     AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
1071     if(!esdevent) {
1072       AliDebug(1, "No esd event available\n");
1073       return;
1074     }
1075
1076     // Check whether primary vertex is available
1077     const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
1078     if(!vtxESDSkip) return;
1079
1080     const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
1081     if(tmptrack){
1082       //case ESD track: take copy constructor
1083
1084       if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1085         vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
1086         isRecalcVertex = kTRUE;
1087       }
1088       if(vtxESDSkip){
1089         AliESDtrack esdtrack(*tmptrack);
1090         fMagField = fEvent->GetMagneticField();
1091         esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD);
1092
1093         if(isRecalcVertex) delete vtxESDSkip;
1094       }
1095     }
1096   }
1097   else {
1098     //case of AOD tracks
1099     AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
1100     if(!aodevent) {
1101       AliDebug(1, "No aod event available\n");
1102       return;
1103     }
1104
1105     // Check whether primary vertex is available
1106     AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
1107     if(!vtxAODSkip) return;
1108
1109     //Case ESD track: take copy constructor
1110     const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
1111     if(tmptrack){
1112
1113       if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1114         vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
1115         isRecalcVertex = kTRUE;
1116       }
1117       if(vtxAODSkip){
1118         AliAODTrack aodtrack(*tmptrack);
1119         AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
1120         fMagField = aodevent->GetMagneticField();
1121         etp.PropagateToDCA(vtxAODSkip, fMagField, kBeampiperadius, dcaD, covD);
1122
1123         if(isRecalcVertex) delete vtxAODSkip;
1124       }
1125     }
1126   }
1127 }
1128
1129 //______________________________________________________
1130 void AliHFEextraCuts::GetHFEImpactParameterCuts(const AliVTrack * const track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
1131         //
1132         // Get HFE impact parameter cut(pt dependent)
1133         //
1134
1135   Double_t pt = track->Pt();
1136   hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt);  // abs R cut
1137   hfeimpactnsigmaRcut = fIPcutParam[3];                                  // sigma cut
1138 }
1139 //______________________________________________________
1140 void AliHFEextraCuts::GetMaxImpactParameterCutR(const AliVTrack * const track, Double_t &maximpactRcut){
1141         //
1142         // Get max impact parameter cut r (pt dependent)
1143         //
1144
1145   Double_t pt = track->Pt();
1146   if(pt > 0.15) {
1147     maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01);  // abs R cut
1148   }
1149   else maximpactRcut = 9999999999.0;
1150 }
1151 //______________________________________________________
1152 void AliHFEextraCuts::GetTOFsignalDxDz(const AliVTrack * const track, Double_t &tofsignalDx, Double_t &tofsignalDz){
1153   //
1154   // TOF matching
1155   //
1156
1157   TString type = track->IsA()->GetName();
1158   if(!type.CompareTo("AliESDtrack")){
1159     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
1160     if(!esdtrack) return;
1161     tofsignalDx = esdtrack->GetTOFsignalDx();
1162     tofsignalDz = esdtrack->GetTOFsignalDz();
1163   }
1164
1165 }
1166
1167 //______________________________________________________
1168 Bool_t AliHFEextraCuts::CheckITSpattern(const AliVTrack *const track) const {
1169   //
1170   // Check if every ITS layer, which has a module which is alive, also
1171   // has an ITS cluster
1172   //
1173   Bool_t patternOK(kTRUE);
1174   Int_t status(0);
1175   for(Int_t ily = 0; ily < 6; ily++){
1176     status = GetITSstatus(track, ily);
1177     if(CheckITSstatus(status)){
1178       // pixel alive, check whether layer has a cluster
1179       if(!TESTBIT(track->GetITSClusterMap(),ily)){
1180         // No cluster even though pixel is alive - reject track
1181         patternOK = kFALSE;
1182         break;
1183       }
1184     }
1185   }
1186   return patternOK;
1187 }
1188
1189 //---------------------------------------------------------------------------
1190 const AliVVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliESDEvent * const esdevent, const AliVTrack * const track) {
1191   //
1192   // This method returns a primary vertex without the daughter tracks of the
1193   // candidate and it recalculates the impact parameters and errors for ESD tracks.
1194   //
1195   // The output vertex is created with "new".
1196   //
1197
1198   //const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
1199
1200   AliVertexerTracks vertexer(fEvent->GetMagneticField());
1201   vertexer.SetITSMode();
1202   vertexer.SetMinClusters(4);
1203   Int_t skipped[2];
1204   skipped[0] = track->GetID();
1205   vertexer.SetSkipTracks(1,skipped);
1206
1207   //diamond constraint
1208   vertexer.SetConstraintOn();
1209   Float_t diamondcovxy[3];
1210   esdevent->GetDiamondCovXY(diamondcovxy);
1211   Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
1212   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1213   AliESDVertex diamond(pos,cov,1.,1);
1214   vertexer.SetVtxStart(&diamond);
1215
1216   const AliVVertex *vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
1217
1218   return vtxESDSkip;
1219 }
1220
1221 //---------------------------------------------------------------------------
1222 AliAODVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliAODEvent * const aod, const AliVTrack * const track) {
1223   //
1224   // This method returns a primary vertex without the daughter tracks of the
1225   // candidate and it recalculates the impact parameters and errors for AOD tracks.
1226   // The output vertex is created with "new".
1227   //
1228
1229   AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
1230   if(!vtxAOD) return 0;
1231   TString title=vtxAOD->GetTitle();
1232   if(!title.Contains("VertexerTracks")) return 0;
1233
1234   AliVertexerTracks vertexer(aod->GetMagneticField());
1235
1236   vertexer.SetITSMode();
1237   vertexer.SetMinClusters(3);
1238   vertexer.SetConstraintOff();
1239
1240   if(title.Contains("WithConstraint")) {
1241     Float_t diamondcovxy[3];
1242     aod->GetDiamondCovXY(diamondcovxy);
1243     Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
1244     Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1245     AliESDVertex diamond(pos,cov,1.,1);
1246     vertexer.SetVtxStart(&diamond);
1247   }
1248   Int_t skipped[2]; for(Int_t i=0;i<2;i++) skipped[i]=-1;
1249   Int_t id = (Int_t)track->GetID();
1250   if(!(id<0)) skipped[0] = id;
1251
1252   /*// exclude tracks with global constrained parameters
1253   Int_t nTracks=aod->GetNumberOfTracks();
1254   for(Int_t i=0; i<nTracks; i++){
1255     t = aod->GetTrack(i);
1256     if(t->TestFilterMask(512)){
1257       id = (Int_t)t->GetID();
1258       skipped[nTrksToSkip++] = id;
1259     }
1260   }*/
1261
1262   vertexer.SetSkipTracks(1,skipped);
1263   AliESDVertex *vtxESDNew = vertexer.FindPrimaryVertex(aod);
1264
1265   if(!vtxESDNew) return 0;
1266   if(vtxESDNew->GetNContributors()<=0) {
1267     delete vtxESDNew; vtxESDNew=NULL;
1268     return 0;
1269   }
1270
1271   // convert to AliAODVertex
1272   Double_t pos[3],cov[6],chi2perNDF;
1273   vtxESDNew->GetXYZ(pos); // position
1274   vtxESDNew->GetCovMatrix(cov); //covariance matrix
1275   chi2perNDF = vtxESDNew->GetChi2toNDF();
1276   delete vtxESDNew; vtxESDNew=NULL;
1277
1278   AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);
1279
1280   return vtxAODNew;
1281 }