PWGGA/GammaConv/AliConversionTrackCuts.cxx

   1 /**************************************************************************
   2  * This file is property of and copyright by the ALICE HLT Project        *
   3  * ALICE Experiment at CERN, All rights reserved.                         *
   4  *                                                                        *
   5  * Primary Author: Svein Lindal <slindal@fys.uio.no>                      *
   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 /// @file   AliConversionTrackCuts.cxx
  17 /// @author Svein Lindal
  18 /// @brief  Base class for analysation of conversion particle - track correlations
  19
  20
  21 #include "AliConversionTrackCuts.h"
  22 //#include "AliAODTrack.h"
  23 #include "AliAODEvent.h"
  24 #include <TFormula.h>
  25 #include <iostream>
  26 #include "TH2F.h"
  27 #include "AliESDtrackCuts.h"
  28
  29 using namespace std;
  30 ClassImp(AliConversionTrackCuts)
  31
  32
  33 const char* AliConversionTrackCuts::fgkCutNames[AliConversionTrackCuts::kNCuts] = {
  34   "nClusTPC",
  35   "FoundFindable",
  36   "Chi2PerNDF",
  37   "Kink",
  38   "DCA_Z",
  39   "DCA_XY",
  40   "TPCRefit"
  41   "kAccTracks"
  42 };
  43
  44
  45
  46 //________________________________________________________________________
  47 AliConversionTrackCuts::AliConversionTrackCuts() :
  48   AliAnalysisCuts(),
  49   fEsdTrackCuts(NULL),
  50   fEsdTrackCutsExtra1(NULL),
  51   fEsdTrackCutsExtra2(NULL),
  52   fEvent(NULL),
  53   fFilterBit(2048),
  54   fDCAZmax(-1),
  55   fDCAXYmax(-1),
  56   fOwnedTracks(),
  57   fInitialized(kFALSE),
  58   fhPhi(NULL),
  59   //  fhPt(NULL),
  60   //fhPhiPt(NULL),
  61   fhdcaxyPt(NULL),
  62   fhdcazPt(NULL),
  63   fhdca(NULL),
  64   fhnclpt(NULL),
  65   fhnclsfpt(NULL),
  66   fhEtaPhi(NULL),
  67   fHistograms(NULL)
  68 {
  69   //Constructor
  70   fOwnedTracks.SetOwner(kTRUE);
  71 }
  72 //________________________________________________________________________
  73 AliConversionTrackCuts::AliConversionTrackCuts(TString name, TString title = "title") :
  74   AliAnalysisCuts(name, title),
  75   fEsdTrackCuts(NULL),
  76   fEsdTrackCutsExtra1(NULL),
  77   fEsdTrackCutsExtra2(NULL),
  78   fEvent(NULL),
  79   fFilterBit(2048),
  80   fDCAZmax(-1),
  81   fDCAXYmax(-1),
  82   fOwnedTracks(),
  83   fInitialized(kFALSE),
  84   fhPhi(NULL),
  85   //fhPt(NULL),
  86   //fhPhiPt(NULL),
  87   fhdcaxyPt(NULL),
  88   fhdcazPt(NULL),
  89   fhdca(NULL),
  90   fhnclpt(NULL),
  91   fhnclsfpt(NULL),
  92   fhEtaPhi(NULL),
  93   fHistograms(NULL)
  94 {
  95   //Constructor
  96   fOwnedTracks.SetOwner(kTRUE);
  97 }
  98
  99
 100 //________________________________________________________________________________
 101  AliConversionTrackCuts::~AliConversionTrackCuts() {
 102    ///destructor
 103    // if(fHistograms)
 104    //    delete fHistograms;
 105    // fHistograms = NULL;
 106
 107    if(fEsdTrackCuts)
 108      delete fEsdTrackCuts;
 109    fEsdTrackCuts = NULL;
 110
 111    if(fEsdTrackCutsExtra1)
 112      delete fEsdTrackCutsExtra1;
 113    fEsdTrackCutsExtra1 = NULL;
 114
 115    if(fEsdTrackCutsExtra2)
 116      delete fEsdTrackCutsExtra2;
 117    fEsdTrackCutsExtra2 = NULL;
 118
 119    fOwnedTracks.Delete();
 120 }
 121
 122 //______________________________________________________________________________
 123 void AliConversionTrackCuts::DefineESDCuts() {
 124   // Reproduces the cuts of the corresponding bit in the ESD->AOD filtering
 125   // (see $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C)
 126   ///Copied from alianalyseleadingue
 127   const Int_t filterbit = fFilterBit;
 128
 129   if (filterbit == 128) {
 130     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 131     fEsdTrackCuts->SetMinNClustersTPC(70);
 132
 133   }  else if (filterbit == 256) {
 134     // syst study
 135     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 136     fEsdTrackCuts->SetMinNClustersTPC(80);
 137     fEsdTrackCuts->SetMaxChi2PerClusterTPC(3);
 138     fEsdTrackCuts->SetMaxDCAToVertexZ(2.7);
 139     fEsdTrackCuts->SetMaxDCAToVertexXY(1.9);
 140
 141   }  else if (filterbit == 512) {
 142     // syst study
 143     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 144     fEsdTrackCuts->SetMinNClustersTPC(60);
 145     fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
 146     fEsdTrackCuts->SetMaxDCAToVertexZ(3.7);
 147     fEsdTrackCuts->SetMaxDCAToVertexXY(2.9);
 148
 149   } else if (filterbit == 1024) {
 150     if(!fEsdTrackCuts) {
 151       fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 152       fEsdTrackCuts->SetMinNClustersTPC(-1);
 153       fEsdTrackCuts->SetMinNCrossedRowsTPC(70);
 154       fEsdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 155     }
 156   } else if (filterbit == 2048)  {
 157     // mimic hybrid tracks
 158     // correspond to esdTrackCutsHTG, but WITHOUT spd constraint. this is checked with the next object
 159     if(!fEsdTrackCuts) {
 160       fEsdTrackCuts = new AliESDtrackCuts();
 161       TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
 162       fEsdTrackCuts->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep, 100);
 163       fEsdTrackCuts->SetMaxChi2PerClusterTPC(4);
 164       fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE);
 165       fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 166       fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
 167       fEsdTrackCuts->SetMaxFractionSharedTPCClusters(0.4);
 168
 169       fEsdTrackCuts->SetMaxDCAToVertexXY(2.4);
 170       fEsdTrackCuts->SetMaxDCAToVertexZ(3.2);
 171       fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
 172
 173       fEsdTrackCuts->SetMaxChi2PerClusterITS(36);
 174       fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 175
 176       fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 177
 178       fEsdTrackCuts->SetEtaRange(-0.9, 0.9);
 179       fEsdTrackCuts->SetPtRange(0.1, 1000000.0);
 180
 181       fEsdTrackCuts->SetRequireITSRefit(kFALSE); //not here, n
 182     }
 183     // Add SPD requirement
 184     fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
 185     fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 186     fEsdTrackCutsExtra1->SetRequireITSRefit(kTRUE);
 187     // A track passing fEsdTrackCuts and fEsdTrackCutsExtra1 corresponds to esdTrackCutsHTG
 188
 189     fEsdTrackCutsExtra2 = new AliESDtrackCuts("No_SPD", "Reject tracks with cluster in SPD");
 190     fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kNone);
 191     // A track passing fEsdTrackCuts and fEsdTrackCutsExtra2 corresponds to esdTrackCutsHTGC and needs to be constrained
 192
 193
 194   }
 195 }
 196
 197 //______________________________________________________________________________
 198 Bool_t AliConversionTrackCuts::AcceptTrack(AliESDtrack * track) {
 199   //Check esd track
 200
 201   if(!fInitialized) {
 202     DefineESDCuts();
 203     if(fDCAXYmax > 0) {
 204       if(fEsdTrackCuts) fEsdTrackCuts->SetMaxDCAToVertexXY(fDCAXYmax);
 205     }
 206     if(fDCAZmax > 0) {
 207       if(fEsdTrackCuts) fEsdTrackCuts->SetMaxDCAToVertexZ(fDCAZmax);
 208     }
 209
 210     fInitialized = kTRUE;
 211   }
 212
 213   FillHistograms(kPreCut, track);
 214
 215
 216
 217   if( !fEsdTrackCuts->IsSelected(track)) return kFALSE;
 218
 219
 220   ///If only one track cuts then it has passed the cuts
 221   if( !(fEsdTrackCutsExtra1 && fEsdTrackCutsExtra2)) {
 222     FillHistograms(1, track);
 223     return kTRUE;
 224   }
 225   ///If passing extra
 226   if (fEsdTrackCutsExtra1 && fEsdTrackCutsExtra1->IsSelected(track)) {
 227     FillHistograms(2, track);
 228     FillHistograms(1, track);
 229     return kTRUE;
 230   }
 231
 232   ///If passing extra2
 233   if (fEsdTrackCutsExtra2 && fEsdTrackCutsExtra2->IsSelected(track)) {
 234     const AliExternalTrackParam * param = track->GetConstrainedParam();
 235     if(param) {
 236       AliESDtrack* esdTrack = new AliESDtrack(*track);
 237       esdTrack->CopyFromVTrack(param);
 238       track = esdTrack;
 239       fOwnedTracks.Add(track);
 240
 241       FillHistograms(3, track);
 242       FillHistograms(1, track);
 243       return kTRUE;
 244     } else {
 245       return kFALSE;
 246     }
 247   } else {
 248     return kFALSE;
 249   }
 250
 251   cout << "error error, should not be herer!"<<endl;
 252   return kFALSE;
 253
 254   // FillHistograms(kPreCut + 1, track);
 255   // return kTRUE;
 256
 257   // fhnclpt->Fill(track->Pt(), track->GetTPCNcls());
 258   // if(track->GetTPCNclsF() > 0) fhnclsfpt->Fill(track->Pt(), ((Double_t) track->GetTPCNcls())/track->GetTPCNclsF());
 259   // FillHistograms(kPreCut + 1, track);
 260
 261   // ///Get impact parameters
 262   // Double_t extCov[15];
 263   // track->GetExternalCovariance(extCov);
 264   // Float_t b[2];
 265   // Float_t bCov[3];
 266   // track->GetImpactParameters(b,bCov);
 267   // if (bCov[0]<=0 || bCov[2]<=0) {
 268   //   AliDebug(1, "Estimated b resolution lower or equal zero!");
 269   //   bCov[0]=0; bCov[2]=0;
 270   // }
 271
 272   // Float_t dcaToVertexXY = b[0];
 273   // Float_t dcaToVertexZ = b[1];
 274   // FillDCAHist(dcaToVertexZ, dcaToVertexXY, track);
 275   // return kTRUE;
 276 }
 277
 278 Bool_t AliConversionTrackCuts::AcceptTrack(AliAODTrack * track) {
 279   //Check aod track
 280
 281   FillHistograms(kPreCut, track);
 282
 283   if(!track->IsHybridGlobalConstrainedGlobal()) return kFALSE;
 284
 285    if (!(track->GetStatus() & AliVTrack::kITSrefit)) {
 286      return kFALSE;
 287    }
 288
 289
 290
 291   //The cluster sharing cut can be done with:
 292   Double_t frac = Double_t(track->GetTPCnclsS()) / Double_t(track->GetTPCncls());
 293   if (frac > 0.4) return kFALSE;
 294
 295
 296   ///Do dca xy cut!
 297   FillHistograms(1, track);
 298
 299   ///DCA
 300   Double_t dca[2] = { -999, -999};
 301   //Bool_t dcaok =
 302   GetDCA(track, dca);
 303   FillDCAHist(dca[1], dca[0], track);
 304
 305
 306   if(track->IsGlobalConstrained()) {
 307     FillHistograms(3, track);
 308   } else {
 309     FillHistograms(2, track);
 310   }
 311
 312   if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
 313
 314   return kTRUE;
 315 }
 316
 317 ///______________________________________________________________________________
 318 Bool_t AliConversionTrackCuts::GetDCA(const AliAODTrack *track, Double_t dca[2]) {
 319   if(track->TestBit(AliAODTrack::kIsDCA)){
 320     dca[0]=track->DCA();
 321     dca[1]=track->ZAtDCA();
 322     return kTRUE;
 323   }
 324
 325   Bool_t ok=kFALSE;
 326   if(fEvent) {
 327     Double_t covdca[3];
 328     //AliAODTrack copy(*track);
 329     AliExternalTrackParam etp; etp.CopyFromVTrack(track);
 330
 331     Float_t xstart = etp.GetX();
 332     if(xstart>3.) {
 333       dca[0]=-999.;
 334     dca[1]=-999.;
 335     //printf("This method can be used only for propagation inside the beam pipe \n");
 336     return kFALSE;
 337     }
 338
 339
 340     AliAODVertex *vtx =(AliAODVertex*)(fEvent->GetPrimaryVertex());
 341     Double_t fBzkG = fEvent->GetMagneticField(); // z componenent of field in kG
 342     ok = etp.PropagateToDCA(vtx,fBzkG,kVeryBig,dca,covdca);
 343     //ok = copy.PropagateToDCA(vtx,fBzkG,kVeryBig,dca,covdca);
 344   }
 345   if(!ok){
 346     dca[0]=-999.;
 347     dca[1]=-999.;
 348   }
 349   return ok;
 350 }
 351
 352
 353
 354 TList * AliConversionTrackCuts::CreateHistograms() {
 355   //Create the histograms
 356
 357   if(!fHistograms) fHistograms = new TList();
 358
 359   fHistograms->SetOwner(kTRUE);
 360   fHistograms->SetName("trackCuts");
 361
 362   fhPhi = new TH2F("phi", "phi", 5, -0.5, 4.5, 32, 0, TMath::TwoPi());
 363   // TAxis * xax = fhPhi->GetXaxis();
 364   // for(Int_t i = 0; i < kNCuts; i++){
 365   //    xax->SetBinLabel(xax->FindFixBin(i), fgkCutNames[i]);
 366   // }
 367   fHistograms->Add(fhPhi);
 368
 369
 370   fhEtaPhi = new TH2F("etahpi", "etaphi", 36, -0.9, 0.9, 32, 0, TMath::TwoPi());
 371   fHistograms->Add(fhEtaPhi);
 372
 373   // fhPt = new TH2F("pt", "pt", kNCuts+2, kPreCut -0.5, kNCuts + 0.5,
 374   //                              20, 0., 20.);
 375   // xax = fhPt->GetXaxis();
 376   // for(Int_t i = 0; i < kNCuts; i++){
 377   //    xax->SetBinLabel(xax->FindFixBin(i), fgkCutNames[i]);
 378   // }
 379   // fHistograms->Add(fhPt);
 380
 381   //  fhPhiPt = new TH2F("phipt", "phipt", 100, 0, 100, 64, 0, TMath::TwoPi());
 382   //fHistograms->Add(fhPhiPt);
 383
 384   fhdcaxyPt = new TH2F("dcaxypt", "dcaxypt", 20, 0, 20, 50, -2.5, 2.5);
 385   fHistograms->Add(fhdcaxyPt);
 386
 387   fhdcazPt = new TH2F("dcazpt", "dcazpt", 20, 0, 20, 70, -3.5, 3.5);
 388   fHistograms->Add(fhdcazPt);
 389
 390   fhdca = new TH2F("dca", "dca", 70, -3.5, 3.5, 50, -2.5, 2.5);
 391   fhdca->SetXTitle("dca z");
 392   fhdca->SetYTitle("dca xy");
 393
 394
 395   fHistograms->Add(fhdca);
 396
 397   // fhnclpt = new TH2F("nclstpcvspt", "nclstpcvspt", 20, 0, 20, 50, 0, 100);
 398   // fHistograms->Add(fhnclpt);
 399
 400   // fhnclsfpt = new TH2F("nclsfpt", "nclsfpt", 20, 0, 20, 60, 0, 1.2);
 401   // fHistograms->Add(fhnclsfpt);
 402
 403   return fHistograms;
 404 }
 405
 406 void AliConversionTrackCuts::FillHistograms(Int_t cutIndex, AliVTrack * track) {
 407
 408   //Fill histograms
 409   if(fhPhi) fhPhi->Fill(cutIndex, track->Phi());
 410   //  if(fhPt) fhPt->Fill(cutIndex, track->Pt());
 411   //if(passed) fhPhiPt->Fill(track->Pt(), track->Phi());
 412
 413 }
 414
 415 void AliConversionTrackCuts::FillDCAHist(Float_t dcaz, Float_t dcaxy, AliVTrack * track) {
 416   if(fhdcaxyPt) fhdcaxyPt->Fill(track->Pt(), dcaxy);
 417   if(fhdcazPt) fhdcazPt->Fill(track->Pt(), dcaz);
 418   if(fhdca) fhdca->Fill(dcaz, dcaxy);
 419 }
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429 //_________________________________________________________________________________________________
 430 void AliConversionTrackCuts::Print(const Option_t *) const
 431 {
 432 //
 433 // Print information on this cut
 434 //
 435
 436   // printf("Cut name                : %s \n", GetName());
 437   // printf("Kink daughters are      : %s \n", (fRejectKinkDaughters ? "rejected" : "accepted"));
 438   // printf("TPC requirements        : clusters/findable %f, min. cluster = %d, max chi2 = %f, %s require refit\n", fTPCClusOverFindable, fTPCminNClusters, fTPCmaxChi2, (fRequireTPCRefit) ? "" : "Don't");
 439   // printf("ITS requirements        : min. cluster = %d (all), %d (SPD), max chi2 = %f \n", fITSminNClusters, fSPDminNClusters, fITSmaxChi2);
 440   // printf("DCA z cut               : fixed to %f cm \n", fDCAZmax);
 441   // printf("DCA xy cut              : fixed to %f cm \n", fDCAXYmax)
 442     ;
 443 }
 444
 445 //_________________________________________________________________________________________________
 446
 447 Bool_t AliConversionTrackCuts::IsSelected(TObject * object ) {
 448   AliAODTrack * aodtrack = dynamic_cast<AliAODTrack*>(object);
 449   if (aodtrack) {
 450     return AcceptTrack(aodtrack);
 451   } else {
 452     AliESDtrack * track = dynamic_cast<AliESDtrack*>(object);
 453     if (track)
 454       return AcceptTrack(track);
 455   }
 456
 457 return kFALSE;
 458 }
 459
 460
 461
 462