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