PWGHF/hfe/AliHFEextraCuts.cxx

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