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 "AliESDEvent.h"
  36 #include "AliESDVertex.h"
  37 #include "AliESDtrack.h"
  38 #include "AliLog.h"
  39 #include "AliMCParticle.h"
  40 #include "AliVEvent.h"
  41 #include "AliVTrack.h"
  42 #include "AliVParticle.h"
  43 #include "AliVertexerTracks.h"
  44 #include "AliVVertex.h"
  45
  46 #include "AliHFEextraCuts.h"
  47
  48 ClassImp(AliHFEextraCuts)
  49
  50 const Int_t AliHFEextraCuts::fgkNQAhistos = 8;
  51
  52 //______________________________________________________
  53 AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
  54   AliCFCutBase(name, title),
  55   fEvent(NULL),
  56   fCutCorrelation(0),
  57   fRequirements(0),
  58   fMinNClustersTPC(0),
  59   fMinNClustersTPCPID(0),
  60   fClusterRatioTPC(0.),
  61   fMinTrackletsTRD(0),
  62   fTRDtrackletsExact(0),
  63   fPixelITS(0),
  64   fTPCclusterDef(0),
  65   fTPCclusterRatioDef(0),
  66   fFractionTPCShared(-1.0),
  67   fCheck(kFALSE),
  68   fQAlist(0x0) ,
  69   fDebugLevel(0)
  70 {
  71   //
  72   // Default Constructor
  73   //
  74   memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
  75   memset(fIPcutParam, 0, sizeof(Float_t) * 4);
  76 }
  77
  78 //______________________________________________________
  79 AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
  80   AliCFCutBase(c),
  81   fEvent(c.fEvent),
  82   fCutCorrelation(c.fCutCorrelation),
  83   fRequirements(c.fRequirements),
  84   fMinNClustersTPC(c.fMinNClustersTPC),
  85   fMinNClustersTPCPID(c.fMinNClustersTPCPID),
  86   fClusterRatioTPC(c.fClusterRatioTPC),
  87   fMinTrackletsTRD(c.fMinTrackletsTRD),
  88   fTRDtrackletsExact(c.fTRDtrackletsExact),
  89   fPixelITS(c.fPixelITS),
  90   fTPCclusterDef(c.fTPCclusterDef),
  91   fTPCclusterRatioDef(c.fTPCclusterRatioDef),
  92   fFractionTPCShared(c.fFractionTPCShared),
  93   fCheck(c.fCheck),
  94   fQAlist(0x0),
  95   fDebugLevel(0)
  96 {
  97   //
  98   // Copy constructor
  99   // Performs a deep copy
 100   //
 101   memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
 102   memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
 103   if(IsQAOn()){
 104     fIsQAOn = kTRUE;
 105     fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
 106     if(fQAlist) fQAlist->SetOwner();
 107   }
 108 }
 109
 110 //______________________________________________________
 111 AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
 112   //
 113   // Assignment operator
 114   //
 115   if(this != &c){
 116     AliCFCutBase::operator=(c);
 117     fEvent = c.fEvent;
 118     fCutCorrelation = c.fCutCorrelation;
 119     fRequirements = c.fRequirements;
 120     fClusterRatioTPC = c.fClusterRatioTPC;
 121     fMinNClustersTPC = c.fMinNClustersTPC;
 122     fMinNClustersTPCPID = c.fMinNClustersTPCPID;
 123     fMinTrackletsTRD = c.fMinTrackletsTRD;
 124     fTRDtrackletsExact = c.fTRDtrackletsExact;
 125     fPixelITS = c.fPixelITS;
 126     fTPCclusterDef = c.fTPCclusterDef;
 127     fTPCclusterRatioDef = c.fTPCclusterRatioDef;
 128     fFractionTPCShared = c.fFractionTPCShared;
 129     fCheck = c.fCheck;
 130     fDebugLevel = c.fDebugLevel;
 131     memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
 132     memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
 133     if(IsQAOn()){
 134       fIsQAOn = kTRUE;
 135       fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
 136       if(fQAlist) fQAlist->SetOwner();
 137     }else fQAlist = 0x0;
 138   }
 139   return *this;
 140 }
 141
 142 //______________________________________________________
 143 AliHFEextraCuts::~AliHFEextraCuts(){
 144   //
 145   // Destructor
 146   //
 147   if(fQAlist) delete fQAlist;
 148 }
 149
 150 //______________________________________________________
 151 void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
 152   //
 153   // Set Virtual event an make a copy
 154   //
 155   if (!event) {
 156     AliError("Pointer to AliVEvent !");
 157     return;
 158   }
 159   TString className(event->ClassName());
 160   if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
 161     AliError("argument must point to an AliESDEvent or AliAODEvent !");
 162     return ;
 163   }
 164   fEvent = (AliVEvent*) event;
 165
 166 }
 167
 168 //______________________________________________________
 169 Bool_t AliHFEextraCuts::IsSelected(TObject *o){
 170   //
 171   // Steering function for the track selection
 172   //
 173   TString type = o->IsA()->GetName();
 174   AliDebug(2, Form("Object type %s", type.Data()));
 175   if(!type.CompareTo("AliESDtrack") || !type.CompareTo("AliAODTrack")){
 176     return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
 177   }
 178   return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
 179 }
 180
 181 //______________________________________________________
 182 Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
 183   //
 184   // Checks cuts on reconstructed tracks
 185   // returns true if track is selected
 186   // QA histograms are filled before track selection and for
 187   // selected tracks after track selection
 188   //
 189   AliDebug(1, "Called");
 190   ULong64_t survivedCut = 0;    // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
 191   if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
 192   // Apply cuts
 193   Float_t impactR, impactZ;
 194   Double_t hfeimpactR, hfeimpactnsigmaR;
 195   Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
 196   Double_t maximpactRcut;
 197   GetImpactParameters(track, impactR, impactZ);
 198   if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
 199     // Protection for PbPb
 200     GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
 201     GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
 202   }
 203   UInt_t nclsTPC = GetTPCncls(track);
 204   UInt_t nclsTPCPID = GetTPCnclusdEdx(track);
 205   // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 206   Float_t fractionSharedClustersTPC = GetTPCsharedClustersRatio(track);
 207   Double_t ratioTPC = GetTPCclusterRatio(track);
 208   UChar_t trdTracklets;
 209   trdTracklets = track->GetTRDntrackletsPID();
 210   UChar_t itsPixel = track->GetITSClusterMap();
 211   Int_t status1 = GetITSstatus(track, 0);
 212   Int_t status2 = GetITSstatus(track, 1);
 213   Bool_t statusL0 = CheckITSstatus(status1);
 214   Bool_t statusL1 = CheckITSstatus(status2);
 215   if(TESTBIT(fRequirements, kTPCfractionShared)) {
 216     // cut on max fraction of shared TPC clusters
 217     if(TMath::Abs(fractionSharedClustersTPC) < fFractionTPCShared) SETBIT(survivedCut, kTPCfractionShared);
 218   }
 219   if(TESTBIT(fRequirements, kMinImpactParamR)){
 220     // cut on min. Impact Parameter in Radial direction
 221     if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
 222   }
 223   if(TESTBIT(fRequirements, kMinImpactParamZ)){
 224     // cut on min. Impact Parameter in Z direction
 225     if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
 226   }
 227   if(TESTBIT(fRequirements, kMaxImpactParamR)){
 228     // cut on max. Impact Parameter in Radial direction
 229     if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
 230   }
 231   if(TESTBIT(fRequirements, kMaxImpactParameterRpar)) {
 232     // HFE Impact parameter cut
 233     GetMaxImpactParameterCutR(track,maximpactRcut);
 234     if(TMath::Abs(impactR) >= maximpactRcut) SETBIT(fRequirements, kMaxImpactParameterRpar);
 235   }
 236   if(TESTBIT(fRequirements, kMaxImpactParamZ)){
 237     // cut on max. Impact Parameter in Z direction
 238     if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
 239   }
 240   if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
 241     // cut on min. HFE Impact Parameter in Radial direction
 242     if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
 243   }
 244   if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
 245     // cut on max. HFE Impact Parameter n sigma in Radial direction
 246     if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
 247   }
 248   if(TESTBIT(fRequirements, kClusterRatioTPC)){
 249     // cut on min ratio of found TPC clusters vs findable TPC clusters
 250     if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
 251   }
 252   if(TESTBIT(fRequirements, kMinTrackletsTRD)){
 253     // cut on minimum number of TRD tracklets
 254     AliDebug(1, Form("Min TRD cut: [%d|%d], Require exact number [%s]\n", fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
 255     if(fTRDtrackletsExact){
 256       if(trdTracklets == fMinTrackletsTRD) {
 257         SETBIT(survivedCut, kMinTrackletsTRD);
 258         AliDebug(1, "Track Selected");
 259       }
 260     }else{
 261       if(trdTracklets >= fMinTrackletsTRD){
 262         SETBIT(survivedCut, kMinTrackletsTRD);
 263         AliDebug(1, "Track Selected");
 264       }
 265       //printf("Min number of tracklets %d\n",fMinTrackletsTRD);
 266     }
 267   }
 268   if(TESTBIT(fRequirements, kMinNClustersTPC)){
 269     // cut on minimum number of TRD tracklets
 270     AliDebug(1, Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
 271     if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
 272   }
 273   if(TESTBIT(fRequirements, kMinNClustersTPCPID)){
 274     AliDebug(1, Form("Min TPC PID cut: [%d|%d]\n", fMinNClustersTPCPID, nclsTPCPID));
 275     if(nclsTPCPID >= fMinNClustersTPCPID) SETBIT(survivedCut, kMinNClustersTPCPID);
 276   }
 277   if(TESTBIT(fRequirements, kPixelITS)){
 278     // cut on ITS pixel layers
 279     AliDebug(1, "ITS cluster Map: ");
 280     //PrintBitMap(itsPixel);
 281     switch(fPixelITS){
 282       case kFirst:
 283         AliDebug(2, "First");
 284               if(TESTBIT(itsPixel, 0))
 285                 SETBIT(survivedCut, kPixelITS);
 286         else
 287                 if(fCheck && !statusL0)
 288                   SETBIT(survivedCut, kPixelITS);
 289                     break;
 290       case kSecond:
 291         AliDebug(2, "Second");
 292               if(TESTBIT(itsPixel, 1))
 293                 SETBIT(survivedCut, kPixelITS);
 294         else
 295                 if(fCheck && !statusL1)
 296                   SETBIT(survivedCut, kPixelITS);
 297                     break;
 298       case kBoth:
 299         AliDebug(2, "Both");
 300               if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
 301                       SETBIT(survivedCut, kPixelITS);
 302               else
 303           if(fCheck && !(statusL0 && statusL1))
 304                         SETBIT(survivedCut, kPixelITS);
 305               break;
 306       case kAny:
 307         AliDebug(2, "Any");
 308               if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
 309                 SETBIT(survivedCut, kPixelITS);
 310         else
 311                 if(fCheck && !(statusL0 || statusL1))
 312                     SETBIT(survivedCut, kPixelITS);
 313                     break;
 314       case kExclusiveSecond:
 315         AliDebug(2, "Exlusive second");
 316         if(fCheck){ // Cut out tracks which pass a dead ITS Layer 0
 317           if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0) && statusL0)
 318             SETBIT(survivedCut, kPixelITS);
 319         } else {
 320           if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0))
 321             SETBIT(survivedCut, kPixelITS);
 322         }
 323         break;
 324       default:
 325         AliDebug(2, "None");
 326         break;
 327     }
 328     AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
 329   }
 330
 331   if(TESTBIT(fRequirements, kTOFPID)){
 332     // cut on TOF pid
 333     if(track->GetStatus() & AliESDtrack::kTOFpid) SETBIT(survivedCut, kTOFPID);
 334   }
 335
 336   if(TESTBIT(fRequirements, kTOFmismatch)){
 337     // cut on TOF mismatch
 338     if(!(track->GetStatus() & AliESDtrack::kTOFmismatch)) SETBIT(survivedCut, kTOFmismatch);
 339   }
 340
 341   if(TESTBIT(fRequirements, kTPCPIDCleanUp)){
 342       // cut on TPC PID cleanup
 343       Int_t fClusterdEdx=GetTPCnclusdEdx(track);
 344       Bool_t fBitsAboveThreshold=GetTPCCountSharedMapBitsAboveThreshold(track);
 345       if((fBitsAboveThreshold==0)&&(fClusterdEdx>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
 346   }
 347
 348   if(TESTBIT(fRequirements, kEMCALmatch)){
 349     if(track->GetStatus() && AliESDtrack::kEMCALmatch) SETBIT(survivedCut, kEMCALmatch);
 350   }
 351
 352   if(fRequirements == survivedCut){
 353     //
 354     // Track selected
 355     //
 356     AliDebug(2, "Track Survived cuts\n");
 357     if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
 358     return kTRUE;
 359   }
 360   AliDebug(2, "Track cut");
 361   if(IsQAOn()) FillCutCorrelation(survivedCut);
 362   return kFALSE;
 363 }
 364
 365 //______________________________________________________
 366 Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
 367   //
 368   // Checks cuts on Monte Carlo tracks
 369   // returns true if track is selected
 370   // QA histograms are filled before track selection and for
 371   // selected tracks after track selection
 372   //
 373   return kTRUE; // not yet implemented
 374 }
 375
 376 //______________________________________________________
 377 void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
 378   //
 379   // Fill the QA histograms for ESD tracks
 380   // Function can be called before cuts or after cut application (second argument)
 381   //
 382   Float_t impactR, impactZ;
 383   GetImpactParameters(track, impactR, impactZ);
 384   TH1 *htmp = NULL;
 385   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * fgkNQAhistos)))) htmp->Fill(impactR);
 386   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * fgkNQAhistos)))) htmp->Fill(impactZ);
 387   // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 388   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * fgkNQAhistos)))) htmp->Fill(GetTPCclusterRatio(track));
 389   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * fgkNQAhistos)))) htmp->Fill(track->GetTRDntrackletsPID());
 390   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * fgkNQAhistos)))) htmp->Fill(GetTPCncls(track));
 391   UChar_t itsPixel = track->GetITSClusterMap();
 392   TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * fgkNQAhistos));
 393   //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
 394   if(pixelHist){
 395     Double_t firstEntry = 0.5;
 396     if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
 397       pixelHist->Fill(firstEntry + 3);
 398     else{
 399       if(itsPixel & BIT(0)){
 400         pixelHist->Fill(firstEntry);
 401         if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
 402         else pixelHist->Fill(firstEntry + 4);
 403       }
 404       if(itsPixel & BIT(1)){
 405         pixelHist->Fill(firstEntry + 1);
 406         if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
 407       }
 408     }
 409   }
 410   // Fill histogram with the status bits
 411   TH1 *hStatusBits = dynamic_cast<TH1 *>(fQAlist->At(6 + when * fgkNQAhistos));
 412   if(hStatusBits) {
 413     hStatusBits->Fill(0);  // Fill first bin with all tracks
 414     if(track->GetStatus() && AliESDtrack::kTOFpid) hStatusBits->Fill(1);
 415     if(!(track->GetStatus() && AliESDtrack::kTOFmismatch)) hStatusBits->Fill(2);
 416     if(track->GetStatus() && AliESDtrack::kEMCALmatch) hStatusBits->Fill(3);
 417     if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
 418   }
 419   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(GetTPCnclusdEdx(track));
 420 }
 421
 422 // //______________________________________________________
 423 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
 424 //   //
 425 //   // Fill the QA histograms for MC tracks
 426 //   // Function can be called before cuts or after cut application (second argument)
 427 //   // Not yet implemented
 428 //   //
 429 // }
 430
 431 //______________________________________________________
 432 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
 433   //
 434   // Fill cut correlation histograms for tracks that didn't pass cuts
 435   //
 436   TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * fgkNQAhistos));
 437   if(correlation){
 438     for(Int_t icut = 0; icut < kNcuts; icut++){
 439       if(!TESTBIT(fRequirements, icut)) continue;
 440       for(Int_t jcut = icut; jcut < kNcuts; jcut++){
 441         if(!TESTBIT(fRequirements, jcut)) continue;
 442         if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
 443                 correlation->Fill(icut, jcut);
 444       }
 445     }
 446   }
 447 }
 448
 449 //______________________________________________________
 450 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
 451   //
 452   // Add QA histograms
 453   // For each cut a histogram before and after track cut is created
 454   // Histos before respectively after cut are stored in different lists
 455   // Additionally a histogram with the cut correlation is created and stored
 456   // in the top directory
 457   //
 458
 459   TH1 *histo1D = 0x0;
 460   TH2 *histo2D = 0x0;
 461   TString cutstr[2] = {"before", "after"};
 462
 463   if(!fQAlist) fQAlist = new TList;  // for internal representation, not owner
 464   for(Int_t icond = 0; icond < 2; icond++){
 465     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * fgkNQAhistos);
 466     fQAlist->AddAt(histo1D, 0 + icond * fgkNQAhistos);
 467     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 468     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 469     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * fgkNQAhistos);
 470     fQAlist->AddAt(histo1D, 1 + icond * fgkNQAhistos);
 471     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 472     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 473     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1)), 2 + icond * fgkNQAhistos);
 474     fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
 475     histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
 476     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 477     qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * fgkNQAhistos);
 478     fQAlist->AddAt(histo1D, 3 + icond * fgkNQAhistos);
 479     histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
 480     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 481     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * fgkNQAhistos);
 482     fQAlist->AddAt(histo1D, 4 + icond * fgkNQAhistos);
 483     histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
 484     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 485     qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * fgkNQAhistos);
 486     fQAlist->AddAt(histo1D, 5 + icond * fgkNQAhistos);
 487     histo1D->GetXaxis()->SetTitle("ITS Pixel");
 488     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 489     Int_t first = histo1D->GetXaxis()->GetFirst();
 490     TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
 491     for(Int_t ilabel = 0; ilabel < 6; ilabel++)
 492       histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
 493     qaList->AddAt((histo1D = new TH1F(Form("%s_trackStatus%s",GetName(),cutstr[icond].Data()), "Track Status", 5, 0, 5)), 6 + icond * fgkNQAhistos);
 494     fQAlist->AddAt(histo1D, 6 + icond * fgkNQAhistos);
 495     histo1D->GetXaxis()->SetTitle("Track Status Bit");
 496     histo1D->GetYaxis()->SetTitle("Number of tracks");
 497     TString tsnames[5] = {"All", "TOFPID", "No TOFmismatch", "EMCALmatch","No TPC shared bit"};
 498     for(Int_t istat = 0; istat < 5; istat++) histo1D->GetXaxis()->SetBinLabel(istat + 1, tsnames[istat].Data());
 499     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);
 500     fQAlist->AddAt(histo1D, 7 + icond * fgkNQAhistos);
 501     histo1D->GetXaxis()->SetTitle("Number of TPC clusters for dEdx calculation");
 502     histo1D->GetYaxis()->SetTitle("counts");
 503   }
 504   // Add cut correlation
 505   qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * fgkNQAhistos);
 506   fQAlist->AddAt(histo2D, 2 * fgkNQAhistos);
 507   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"};
 508   Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
 509   for(Int_t icut = 0; icut < kNcuts; icut++){
 510     histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
 511     histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
 512   }
 513 }
 514
 515 //______________________________________________________
 516 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
 517   for(Int_t ibit = 32; ibit--; )
 518     printf("%d", bitmap & BIT(ibit) ? 1 : 0);
 519   printf("\n");
 520 }
 521
 522 //______________________________________________________
 523 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
 524   //
 525   // Check whether ITS area is dead
 526   //
 527   Bool_t status;
 528   switch(itsStatus){
 529     case 2: status = kFALSE; break;
 530     case 3: status = kFALSE; break;
 531     case 7: status = kFALSE; break;
 532     default: status = kTRUE;
 533   }
 534   return status;
 535 }
 536
 537 //______________________________________________________
 538 Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
 539         //
 540         // Check ITS layer status
 541         //
 542         Int_t status = 0;
 543         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 544                 Int_t det;
 545                 Float_t xloc, zloc;
 546                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 547                 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
 548         }
 549         return status;
 550 }
 551
 552
 553 //______________________________________________________
 554 Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track){
 555   //
 556   // Checks if number of shared bits is above 1
 557   //
 558   Int_t nsharebit = 1;
 559   TString type = track->IsA()->GetName();
 560   if(!type.CompareTo("AliESDtrack")){
 561     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 562     if(esdtrack){ // coverity
 563         TBits shared = esdtrack->GetTPCSharedMap();
 564         if(shared.CountBits() >= 2) nsharebit=1;
 565         else nsharebit=0;
 566     }
 567   }
 568   return nsharebit;
 569 }
 570
 571
 572
 573 //______________________________________________________
 574 UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
 575         //
 576         // Get Number of findable clusters in the TPC
 577         //
 578         Int_t nClusters = 0; // in case no Information available consider all clusters findable
 579         TString type = track->IsA()->GetName();
 580         if(!type.CompareTo("AliESDtrack")){
 581     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 582     if(esdtrack){ // coverity
 583       if(TESTBIT(fTPCclusterDef, kFoundIter1)){
 584                     nClusters = esdtrack->GetTPCNclsIter1();
 585         AliDebug(2, ("Using def kFoundIter1"));
 586       } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
 587         AliDebug(2, ("Using def kCrossedRows"));
 588         nClusters = static_cast<UInt_t>(esdtrack->GetTPCClusterInfo(2,1));
 589       } else{
 590         AliDebug(2, ("Using def kFound"));
 591                     nClusters = esdtrack->GetTPCNcls();
 592       }
 593     }
 594   }
 595         else if(!type.CompareTo("AliAODTrack")){
 596                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 597     if(aodtrack){
 598             const TBits &tpcmap = aodtrack->GetTPCClusterMap();
 599             for(UInt_t ibit = 0; ibit < tpcmap.GetNbits(); ibit++)
 600                     if(tpcmap.TestBitNumber(ibit)) nClusters++;
 601     }
 602         }
 603         return nClusters;
 604 }
 605
 606 //______________________________________________________
 607 UInt_t AliHFEextraCuts::GetTPCnclusdEdx(AliVTrack *track){
 608   //
 609   // Get number of TPC cluster used to calculate dEdx
 610   //
 611   Int_t nClustersdEdx = 0;
 612   TString type = track->IsA()->GetName();
 613   if(!type.CompareTo("AliESDtrack")){
 614     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 615     if(esdtrack){ // coverity
 616       nClustersdEdx = esdtrack->GetTPCsignalN();
 617     }
 618   }
 619   return nClustersdEdx;
 620 }
 621
 622
 623
 624 //______________________________________________________
 625 Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
 626   //
 627   // Get fraction of shared TPC clusters
 628   //
 629   Float_t fracClustersTPCShared = 0.0;
 630   TString type = track->IsA()->GetName();
 631   if(!type.CompareTo("AliESDtrack")){
 632     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 633     if(esdtrack){ // coverity
 634       Float_t nClustersTPC = esdtrack->GetTPCclusters(0);
 635       Int_t nClustersTPCShared = esdtrack->GetTPCnclsS();
 636       if (nClustersTPC!=0) {
 637                fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
 638       }
 639     }
 640   }
 641   return fracClustersTPCShared;
 642 }
 643
 644 //______________________________________________________
 645 Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
 646   //
 647   // Get Ratio of found / findable clusters for different definitions
 648   // Only implemented for ESD tracks
 649   //
 650   Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
 651   TString type = track->IsA()->GetName();
 652   if(!type.CompareTo("AliESDtrack")){
 653     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 654     if(esdtrack){ // coverity
 655       if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
 656         AliDebug(2, "Using ratio def kCROverFindable");
 657         clusterRatio = esdtrack->GetTPCNclsF() ? esdtrack->GetTPCClusterInfo(2,1)/static_cast<Double_t>(esdtrack->GetTPCNclsF()) : 1.; // crossed rows/findable
 658       } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
 659         AliDebug(2, "Using ratio def kFoundOverCR");
 660         clusterRatio = esdtrack->GetTPCClusterInfo(2,0); // found/crossed rows
 661       } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
 662         AliDebug(2, "Using ratio def kFoundOverFindableIter1");
 663         clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
 664       } else {
 665         AliDebug(2, "Using ratio def kFoundOverFindable");
 666         clusterRatio = esdtrack->GetTPCNclsF() ? static_cast<Double_t>(esdtrack->GetTPCNcls())/static_cast<Double_t>(esdtrack->GetTPCNclsF()) : 1.; // found/findable
 667       }
 668     }
 669   }
 670   return clusterRatio;
 671 }
 672
 673 //______________________________________________________
 674 void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
 675         //
 676         // Get impact parameter
 677         //
 678         TString type = track->IsA()->GetName();
 679         if(!type.CompareTo("AliESDtrack")){
 680                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 681     if(esdtrack) esdtrack->GetImpactParameters(radial, z);
 682         }
 683         else if(!type.CompareTo("AliAODTrack")){
 684                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 685     if(aodtrack){
 686                   Double_t xyz[3];
 687                   aodtrack->XYZAtDCA(xyz);
 688                   z = xyz[2];
 689                   radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
 690     }
 691         }
 692 }
 693
 694 //______________________________________________________
 695 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
 696         //
 697         // Get HFE impact parameter (with recalculated primary vertex)
 698         //
 699         dcaxy=0;
 700         dcansigmaxy=0;
 701   if(!fEvent){
 702     AliDebug(1, "No Input event available\n");
 703     return;
 704   }
 705   const Double_t kBeampiperadius=3.;
 706   TString type = track->IsA()->GetName();
 707   Double_t dcaD[2]={-999.,-999.};
 708   Double_t covD[3]={-999.,-999.,-999.};
 709   Float_t dcaF[2]={-999.,-999.};
 710   Float_t covF[3]={-999.,-999.,-999.};
 711
 712   AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
 713   if(!esdevent) {
 714     AliDebug(1, "No esd event available\n");
 715     return;
 716   }
 717   const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
 718   if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
 719    // recalculate primary vertex for peri. and pp
 720    AliVertexerTracks vertexer(fEvent->GetMagneticField());
 721    vertexer.SetITSMode();
 722    vertexer.SetMinClusters(4);
 723          Int_t skipped[2];
 724    skipped[0] = track->GetID();
 725    vertexer.SetSkipTracks(1,skipped);
 726
 727  //----diamond constraint-----------------------------
 728    vertexer.SetConstraintOn();
 729    Float_t diamondcovxy[3];
 730    esdevent->GetDiamondCovXY(diamondcovxy);
 731    Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
 732    Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 733    AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
 734    vertexer.SetVtxStart(diamond);
 735    delete diamond; diamond=NULL;
 736  //----------------------------------------------------
 737
 738    vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
 739
 740    // Getting the DCA
 741    // Propagation always done on a working copy to not disturb the track params of the original track
 742    AliESDtrack *esdtrack = NULL;
 743    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 744      // Case ESD track: take copy constructor
 745      AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 746      if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 747    } else {
 748     // Case AOD track: take different constructor
 749     esdtrack = new AliESDtrack(track);
 750    }
 751    if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD)){
 752     // protection
 753     dcaxy = dcaD[0];
 754     if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
 755     if(!covD[0]) dcansigmaxy = -49.;
 756    }
 757    delete esdtrack;
 758    delete vtxESDSkip;
 759   }
 760   else{
 761    AliESDtrack *esdtrack = NULL;
 762    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 763     // Case ESD track: take copy constructor
 764     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 765     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 766    } else {
 767     // Case AOD track: take different constructor
 768     esdtrack = new AliESDtrack(track);
 769    }
 770    if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF);
 771    dcaxy = dcaF[0];
 772    if(covF[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covF[0]);
 773    if(!covF[0]) dcansigmaxy = -49.;
 774    delete esdtrack;
 775   }
 776 }
 777
 778
 779 //______________________________________________________
 780 void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
 781         //
 782         // Get HFE impact parameter cut(pt dependent)
 783         //
 784
 785   TString type = track->IsA()->GetName();
 786   if(!type.CompareTo("AliESDtrack")){
 787     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 788     if(!esdtrack) return;
 789     Double_t pt = esdtrack->Pt();
 790     hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt);  // abs R cut
 791     hfeimpactnsigmaRcut = fIPcutParam[3];                                  // sigma cut
 792   }
 793 }
 794 //______________________________________________________
 795 void AliHFEextraCuts::GetMaxImpactParameterCutR(AliVTrack *track, Double_t &maximpactRcut){
 796         //
 797         // Get max impact parameter cut r (pt dependent)
 798         //
 799
 800   TString type = track->IsA()->GetName();
 801   if(!type.CompareTo("AliESDtrack")){
 802     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 803     if(!esdtrack) return;
 804     Double_t pt = esdtrack->Pt();
 805     if(pt > 0.15) {
 806       maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01);  // abs R cut
 807     }
 808     else maximpactRcut = 9999999999.0;
 809   }
 810 }