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 - track cuts
  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(3.2*3.2),
  55   fDCAXYmax(2.4*2.4),
  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     if(!fEsdTrackCuts) {
 131       fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 132       fEsdTrackCuts->SetMinNClustersTPC(70);
 133     }
 134   }  else if (filterbit == 256) {
 135     if(!fEsdTrackCuts) {
 136       // syst study
 137       fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 138       fEsdTrackCuts->SetMinNClustersTPC(80);
 139       fEsdTrackCuts->SetMaxChi2PerClusterTPC(3);
 140       fEsdTrackCuts->SetMaxDCAToVertexZ(2.7);
 141       fEsdTrackCuts->SetMaxDCAToVertexXY(1.9);
 142     }
 143   }  else if (filterbit == 512) {
 144     if(!fEsdTrackCuts) {
 145       // syst study
 146       fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 147       fEsdTrackCuts->SetMinNClustersTPC(60);
 148       fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
 149       fEsdTrackCuts->SetMaxDCAToVertexZ(3.7);
 150       fEsdTrackCuts->SetMaxDCAToVertexXY(2.9);
 151     }
 152   } else if (filterbit == 1024) {
 153     if(!fEsdTrackCuts) {
 154       fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 155       fEsdTrackCuts->SetMinNClustersTPC(-1);
 156       fEsdTrackCuts->SetMinNCrossedRowsTPC(70);
 157       fEsdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 158     }
 159   } else if (filterbit == 2048)  {
 160     // mimic hybrid tracks
 161     // correspond to esdTrackCutsHTG, but WITHOUT spd constraint. this is checked with the next object
 162     if(!fEsdTrackCuts) {
 163       fEsdTrackCuts = new AliESDtrackCuts();
 164       TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
 165       fEsdTrackCuts->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep, 100);
 166       fEsdTrackCuts->SetMaxChi2PerClusterTPC(4);
 167       fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE);
 168       fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 169       fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
 170       fEsdTrackCuts->SetMaxFractionSharedTPCClusters(0.4);
 171
 172       fEsdTrackCuts->SetMaxDCAToVertexXY(2.4);
 173       fEsdTrackCuts->SetMaxDCAToVertexZ(3.2);
 174       fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
 175
 176       fEsdTrackCuts->SetMaxChi2PerClusterITS(36);
 177       fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 178
 179       fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 180
 181       fEsdTrackCuts->SetEtaRange(-0.9, 0.9);
 182       fEsdTrackCuts->SetPtRange(0.1, 1000000.0);
 183
 184       fEsdTrackCuts->SetRequireITSRefit(kFALSE); //not here, n
 185     }
 186     // Add SPD requirement
 187     fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
 188     fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 189     fEsdTrackCutsExtra1->SetRequireITSRefit(kTRUE);
 190     // A track passing fEsdTrackCuts and fEsdTrackCutsExtra1 corresponds to esdTrackCutsHTG
 191
 192     fEsdTrackCutsExtra2 = new AliESDtrackCuts("No_SPD", "Reject tracks with cluster in SPD");
 193     fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kNone);
 194     // A track passing fEsdTrackCuts and fEsdTrackCutsExtra2 corresponds to esdTrackCutsHTGC and needs to be constrained
 195
 196
 197   }
 198 }
 199
 200
 201 //______________________________________________________________________________
 202 Bool_t AliConversionTrackCuts::AcceptTrack(AliESDtrack * track) {
 203   //Check esd track
 204   FillHistograms(kPreCut, track);
 205
 206
 207   if( fFilterBit == 256) {
 208
 209     ///Standalone tpc tracks constrained
 210     const AliExternalTrackParam * param = track->GetConstrainedParam();
 211     if(param) {
 212       AliESDtrack* esdTrack = new AliESDtrack(*track);
 213       esdTrack->CopyFromVTrack(param);
 214       track = esdTrack;
 215       fOwnedTracks.Add(track);
 216
 217       if( !fEsdTrackCuts->IsSelected(track)) return kFALSE;
 218
 219       FillHistograms(1, track);
 220
 221       Double_t dca[2];
 222       GetDCA(track, dca);
 223
 224       FillDCAHist(dca[1], dca[0], track);
 225       if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
 226
 227       return kTRUE;
 228     } else {
 229       return kFALSE;
 230     }
 231
 232     return kFALSE;
 233   }
 234
 235
 236   if(!fInitialized) {
 237     DefineESDCuts();
 238     // if(fDCAXYmax > 0) {
 239     //   if(fEsdTrackCuts) fEsdTrackCuts->SetMaxDCAToVertexXY(fDCAXYmax);
 240     // }
 241     // if(fDCAZmax > 0) {
 242     //   if(fEsdTrackCuts) fEsdTrackCuts->SetMaxDCAToVertexZ(fDCAZmax);
 243     // }
 244
 245     fInitialized = kTRUE;
 246   }
 247
 248
 249   Double_t dca[2];
 250   GetDCA(track, dca);
 251
 252
 253   ///If only one track cuts then it has passed the cuts
 254   if( !(fEsdTrackCutsExtra1 && fEsdTrackCutsExtra2)) {
 255     FillHistograms(1, track);
 256     FillDCAHist(dca[1], dca[0], track);
 257     if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
 258     return kTRUE;
 259   }
 260
 261   ///If passing extra
 262   if (fEsdTrackCutsExtra1 && fEsdTrackCutsExtra1->IsSelected(track)) {
 263     FillHistograms(1, track);
 264     FillHistograms(2, track);
 265
 266     FillDCAHist(dca[1], dca[0], track);
 267     if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
 268
 269     return kTRUE;
 270   }
 271
 272   ///If passing extra2
 273   if (fEsdTrackCutsExtra2 && fEsdTrackCutsExtra2->IsSelected(track)) {
 274     const AliExternalTrackParam * param = track->GetConstrainedParam();
 275     if(param) {
 276       AliESDtrack* esdTrack = new AliESDtrack(*track);
 277       esdTrack->CopyFromVTrack(param);
 278       track = esdTrack;
 279       fOwnedTracks.Add(track);
 280
 281       FillHistograms(3, track);
 282       FillHistograms(1, track);
 283
 284       FillDCAHist(dca[1], dca[0], track);
 285       if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
 286
 287       return kTRUE;
 288     } else {
 289       return kFALSE;
 290     }
 291   } else {
 292     return kFALSE;
 293   }
 294
 295   cout << "error error, should not be herer!"<<endl;
 296   return kFALSE;
 297
 298   // FillHistograms(kPreCut + 1, track);
 299   // return kTRUE;
 300
 301   // fhnclpt->Fill(track->Pt(), track->GetTPCNcls());
 302   // if(track->GetTPCNclsF() > 0) fhnclsfpt->Fill(track->Pt(), ((Double_t) track->GetTPCNcls())/track->GetTPCNclsF());
 303   // FillHistograms(kPreCut + 1, track);
 304
 305   // ///Get impact parameters
 306   // Double_t extCov[15];
 307   // track->GetExternalCovariance(extCov);
 308   // return kTRUE;
 309 }
 310
 311 Bool_t AliConversionTrackCuts::AcceptTrack(AliAODTrack * track) {
 312   //Check aod track
 313
 314   FillHistograms(kPreCut, track);
 315
 316   if (fFilterBit == 768) {
 317     if(!track->IsHybridGlobalConstrainedGlobal()) return kFALSE;
 318
 319     if (!(track->GetStatus() & AliVTrack::kITSrefit)) {
 320       return kFALSE;
 321     }
 322
 323     //The cluster sharing cut can be done with:
 324     Double_t frac = Double_t(track->GetTPCnclsS()) / Double_t(track->GetTPCncls());
 325     if (frac > 0.4) return kFALSE;
 326
 327     ///Do dca xy cut!
 328     FillHistograms(1, track);
 329
 330     ///DCA
 331     Double_t dca[2] = { -999, -999};
 332     //Bool_t dcaok =
 333     GetDCA(track, dca);
 334     FillDCAHist(dca[1], dca[0], track);
 335
 336
 337     if(track->IsGlobalConstrained()) {
 338       FillHistograms(3, track);
 339     } else {
 340       FillHistograms(2, track);
 341     }
 342
 343     if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
 344
 345     return kTRUE;
 346
 347     ////////////////////////////////
 348     //// Standalone
 349     ////////////////////////////////
 350   } else  if(fFilterBit == 256) {
 351     if(!track->IsTPCConstrained()) return kFALSE;
 352
 353
 354
 355     ///DCA
 356     Double_t dca[2] = { -999, -999};
 357     GetDCA(track, dca);
 358
 359     if( (dca[0]*dca[0]/fDCAXYmax + dca[1]*dca[1]/fDCAZmax) > 1 ) {
 360       FillHistograms(3, track);
 361       return kFALSE;
 362     }
 363
 364     if(track->GetTPCncls() < 70) {
 365       FillHistograms(4, track);
 366       return kFALSE;
 367     }
 368
 369     AliAODVertex * vtx = track->GetProdVertex();
 370     if (vtx->GetType() == AliAODVertex::kKink ) {
 371       FillHistograms(5, track);
 372       return kFALSE;
 373     }
 374
 375     if(track->Chi2perNDF() > 36) {
 376       FillHistograms(6, track);
 377       return kFALSE;
 378     }
 379     if(track->Chi2perNDF() > 26) {
 380       FillHistograms(7, track);
 381       return kFALSE;
 382     }
 383     if(track->Chi2perNDF() > 16) {
 384       FillHistograms(8, track);
 385       return kFALSE;
 386     }
 387     if(track->Chi2perNDF() > 4) {
 388       FillHistograms(9, track);
 389       return kFALSE;
 390     }
 391
 392
 393
 394     FillDCAHist(dca[1], dca[0], track);
 395
 396     FillHistograms(2, track);
 397     if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
 398     return kTRUE;
 399
 400   }
 401   return kFALSE;
 402 }
 403
 404
 405 ///______________________________________________________________________________
 406 Bool_t AliConversionTrackCuts::GetDCA(const AliESDtrack *track, Double_t dcaxyz[2]) {
 407   ///Get track dca esd trck
 408   Float_t dca[2];
 409   Float_t bCov[3];
 410   track->GetImpactParameters(dca,bCov);
 411   if (bCov[0]<=0 || bCov[2]<=0) {
 412     AliDebug(1, "Estimated b resolution lower or equal zero!");
 413     bCov[0]=0; bCov[2]=0;
 414     return kFALSE;
 415   }
 416
 417   dcaxyz[0] = dca[0];
 418   dcaxyz[1] = dca[1];
 419
 420   return kTRUE;
 421 }
 422
 423 ///_____________________________________________________________________________
 424 Bool_t AliConversionTrackCuts::GetDCA(const AliAODTrack *track, Double_t dca[2]) {
 425   ///Get track dca aod trck
 426   if(track->TestBit(AliAODTrack::kIsDCA)){
 427     dca[0]=track->DCA();
 428     dca[1]=track->ZAtDCA();
 429     return kTRUE;
 430   }
 431
 432   Bool_t ok=kFALSE;
 433   if(fEvent) {
 434     Double_t covdca[3];
 435     //AliAODTrack copy(*track);
 436     AliExternalTrackParam etp; etp.CopyFromVTrack(track);
 437
 438     Float_t xstart = etp.GetX();
 439     if(xstart>3.) {
 440       dca[0]=-999.;
 441       dca[1]=-999.;
 442     //printf("This method can be used only for propagation inside the beam pipe \n");
 443     return kFALSE;
 444     }
 445
 446
 447     AliAODVertex *vtx =(AliAODVertex*)(fEvent->GetPrimaryVertex());
 448     Double_t fBzkG = fEvent->GetMagneticField(); // z componenent of field in kG
 449     ok = etp.PropagateToDCA(vtx,fBzkG,kVeryBig,dca,covdca);
 450     //ok = copy.PropagateToDCA(vtx,fBzkG,kVeryBig,dca,covdca);
 451   }
 452   if(!ok){
 453     dca[0]=-999.;
 454     dca[1]=-999.;
 455   }
 456   return ok;
 457 }
 458
 459
 460
 461 TList * AliConversionTrackCuts::CreateHistograms() {
 462   //Create the histograms
 463
 464   if(!fHistograms) fHistograms = new TList();
 465
 466   fHistograms->SetOwner(kTRUE);
 467   fHistograms->SetName("trackCuts");
 468
 469   fhPhi = new TH2F(Form("phi_%s", GetName()), Form("phi_%s", GetTitle()), 5, -0.5, 4.5, 32, 0, TMath::TwoPi());
 470   // TAxis * xax = fhPhi->GetXaxis();
 471   // for(Int_t i = 0; i < kNCuts; i++){
 472   //    xax->SetBinLabel(xax->FindFixBin(i), fgkCutNames[i]);
 473   // }
 474   fHistograms->Add(fhPhi);
 475
 476
 477   fhEtaPhi = new TH2F(Form("etaphi_%s",GetName()), Form("etaphi_%s", GetTitle()), 36, -0.9, 0.9, 32, 0, TMath::TwoPi());
 478   fHistograms->Add(fhEtaPhi);
 479
 480   // fhPt = new TH2F("pt", "pt", kNCuts+2, kPreCut -0.5, kNCuts + 0.5,
 481   //                              20, 0., 20.);
 482   // xax = fhPt->GetXaxis();
 483   // for(Int_t i = 0; i < kNCuts; i++){
 484   //    xax->SetBinLabel(xax->FindFixBin(i), fgkCutNames[i]);
 485   // }
 486   // fHistograms->Add(fhPt);
 487
 488   //  fhPhiPt = new TH2F("phipt", "phipt", 100, 0, 100, 64, 0, TMath::TwoPi());
 489   //fHistograms->Add(fhPhiPt);
 490
 491   fhdcaxyPt = new TH2F(Form("dcaxypt_%s", GetName()),  Form("dcaxypt_%s", GetTitle()), 20, 0, 20, 50, -2.5, 2.5);
 492   fHistograms->Add(fhdcaxyPt);
 493
 494   fhdcazPt = new TH2F(Form("dcazpt_%s", GetName()),  Form("dcazpt_%s", GetTitle()), 20, 0, 20, 70, -3.5, 3.5);
 495   fHistograms->Add(fhdcazPt);
 496
 497   fhdca = new TH2F(Form("dca_%s", GetName()),  Form("dca_%s", GetTitle()), 70, -3.5, 3.5, 50, -2.5, 2.5);
 498   fhdca->SetXTitle("dca z");
 499   fhdca->SetYTitle("dca xy");
 500
 501
 502   fHistograms->Add(fhdca);
 503
 504   // fhnclpt = new TH2F("nclstpcvspt", "nclstpcvspt", 20, 0, 20, 50, 0, 100);
 505   // fHistograms->Add(fhnclpt);
 506
 507   // fhnclsfpt = new TH2F("nclsfpt", "nclsfpt", 20, 0, 20, 60, 0, 1.2);
 508   // fHistograms->Add(fhnclsfpt);
 509
 510   return fHistograms;
 511 }
 512
 513 void AliConversionTrackCuts::FillHistograms(Int_t cutIndex, AliVTrack * track) {
 514
 515   //Fill histograms
 516   if(fhPhi) fhPhi->Fill(cutIndex, track->Phi());
 517   //  if(fhPt) fhPt->Fill(cutIndex, track->Pt());
 518   //if(passed) fhPhiPt->Fill(track->Pt(), track->Phi());
 519
 520 }
 521
 522 void AliConversionTrackCuts::FillDCAHist(Float_t dcaz, Float_t dcaxy, AliVTrack * track) {
 523   if(fhdcaxyPt) fhdcaxyPt->Fill(track->Pt(), dcaxy);
 524   if(fhdcazPt) fhdcazPt->Fill(track->Pt(), dcaz);
 525   if(fhdca) fhdca->Fill(dcaz, dcaxy);
 526 }
 527
 528
 529
 530
 531 //_________________________________________________________________________________________________
 532 void AliConversionTrackCuts::Print(const Option_t *) const
 533 {
 534 //
 535 // Print information on this cut
 536 //
 537
 538   // printf("Cut name                : %s \n", GetName());
 539   // printf("Kink daughters are      : %s \n", (fRejectKinkDaughters ? "rejected" : "accepted"));
 540   // printf("TPC requirements        : clusters/findable %f, min. cluster = %d, max chi2 = %f, %s require refit\n", fTPCClusOverFindable, fTPCminNClusters, fTPCmaxChi2, (fRequireTPCRefit) ? "" : "Don't");
 541   // printf("ITS requirements        : min. cluster = %d (all), %d (SPD), max chi2 = %f \n", fITSminNClusters, fSPDminNClusters, fITSmaxChi2);
 542   // printf("DCA z cut               : fixed to %f cm \n", fDCAZmax);
 543   // printf("DCA xy cut              : fixed to %f cm \n", fDCAXYmax)
 544     ;
 545 }
 546
 547 //_________________________________________________________________________________________________
 548
 549 Bool_t AliConversionTrackCuts::IsSelected(TObject * object ) {
 550   AliAODTrack * aodtrack = dynamic_cast<AliAODTrack*>(object);
 551   if (aodtrack) {
 552     return AcceptTrack(aodtrack);
 553   } else {
 554     AliESDtrack * track = dynamic_cast<AliESDtrack*>(object);
 555     if (track)
 556       return AcceptTrack(track);
 557   }
 558
 559 return kFALSE;
 560 }
 561
 562
 563
 564