PWGHF/correlationHF/AliDxHFEParticleSelectionD0.cxx

   1 // $Id$
   2
   3 //**************************************************************************
   4 //* This file is property of and copyright by the ALICE Project            *
   5 //* ALICE Experiment at CERN, All rights reserved.                         *
   6 //*                                                                        *
   7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
   8 //*                  Sedat Altinpinar <Sedat.Altinpinar@cern.ch>           *
   9 //*                  Hege Erdal       <hege.erdal@gmail.com>               *
  10 //*                                                                        *
  11 //* Permission to use, copy, modify and distribute this software and its   *
  12 //* documentation strictly for non-commercial purposes is hereby granted   *
  13 //* without fee, provided that the above copyright notice appears in all   *
  14 //* copies and that both the copyright notice and this permission notice   *
  15 //* appear in the supporting documentation. The authors make no claims     *
  16 //* about the suitability of this software for any purpose. It is          *
  17 //* provided "as is" without express or implied warranty.                  *
  18 //**************************************************************************
  19
  20 /// @file   AliDxHFEParticleSelectionD0.cxx
  21 /// @author Sedat Altinpinar, Hege Erdal, Matthias Richter
  22 /// @date   2012-03-19
  23 /// @brief  D0 selection for D0-HFE correlation
  24 ///
  25
  26 #include "AliDxHFEParticleSelectionD0.h"
  27 #include "AliVParticle.h"
  28 //#include "AliAnalysisCuts.h"       // required dependency libANALYSISalice.so
  29 //#include "AliFlowTrackSimple.h"    // required dependency libPWGflowBase.so
  30 //#include "AliFlowCandidateTrack.h" // required dependency libPWGflowTasks.so
  31 //#include "AliCFContainer.h"        // required dependency libCORRFW.so
  32 #include "AliAODRecoDecayHF2Prong.h" // libPWGHFvertexingHF
  33 #include "AliRDHFCutsD0toKpi.h"
  34 #include "TObjArray.h"
  35 #include "THnSparse.h"
  36 #include "AliReducedParticle.h"
  37 #include "TAxis.h"
  38 #include "TString.h"
  39 #include <iostream>
  40 #include <cerrno>
  41 #include <memory>
  42
  43 using namespace std;
  44
  45 /// ROOT macro for the implementation of ROOT specific class methods
  46 ClassImp(AliDxHFEParticleSelectionD0)
  47
  48 AliDxHFEParticleSelectionD0::AliDxHFEParticleSelectionD0(const char* opt)
  49   : AliDxHFEParticleSelection("D0", opt)
  50   , fD0Properties(NULL)
  51   , fD0Daughter0(NULL)
  52   , fD0Daughter1(NULL)
  53   , fCuts(NULL)
  54   , fFillOnlyD0D0bar(0)
  55   , fD0InvMass(0.0)
  56   , fPtBin(-1)
  57   , fHistoList(NULL)
  58 {
  59   // constructor
  60   //
  61   //
  62   //
  63   //
  64   ParseArguments(opt);
  65 }
  66
  67 AliDxHFEParticleSelectionD0::~AliDxHFEParticleSelectionD0()
  68 {
  69   // destructor
  70   if (fD0Properties) {
  71     delete fD0Properties;
  72     fD0Properties=NULL;
  73   }
  74   if (fD0Daughter0) {
  75     delete fD0Daughter0;
  76     fD0Daughter0=NULL;
  77   }
  78   if (fD0Daughter1) {
  79     delete fD0Daughter1;
  80     fD0Daughter1=NULL;
  81   }
  82   if (fHistoList){
  83     delete fHistoList;
  84     fHistoList=NULL;
  85   }
  86
  87   // Note: external object deleted elsewhere
  88   fCuts=NULL;
  89 }
  90
  91 const char* AliDxHFEParticleSelectionD0::fgkDgTrackControlBinNames[]={
  92   "Pt",
  93   "Phi",
  94   "Ptbin",
  95   "D0InvMass",
  96   "Eta"
  97 };
  98
  99 const char* AliDxHFEParticleSelectionD0::fgkCutBinNames[]={
 100   "nDstar->D0",
 101   "nCandSel(Tr)",
 102   "IsInFiducialAcceptance",
 103   "ptbin-1",
 104   "No daugthers",
 105   "Selectioncode 0",
 106   "Selected D0",
 107   "Selected D0bar",
 108   "Selected as both"
 109 };
 110
 111
 112 int AliDxHFEParticleSelectionD0::InitControlObjects()
 113 {
 114   /// init the control objects, can be overloaded by childs which should
 115   /// call AliDxHFEParticleSelection::InitControlObjects() explicitly
 116
 117   fD0Properties=DefineTHnSparse();
 118   AddControlObject(fD0Properties);
 119
 120   //Adding control objects for the daughters
 121   InitControlObjectsDaughters("pi information",0);
 122   InitControlObjectsDaughters("K information",1);
 123   AliInfo(Form("D0 filling scheme: %d\n",fFillOnlyD0D0bar));
 124
 125   fHistoList=new TList;
 126   fHistoList->SetName("D0 Histograms");
 127   fHistoList->SetOwner();
 128
 129   // Histogram storing which cuts have been applied to the tracks
 130   fHistoList->Add(CreateControlHistogram("fWhichCutD0","effective cut for a rejected particle", kNCutLabels, fgkCutBinNames));
 131
 132   AddControlObject(fHistoList);
 133
 134   return AliDxHFEParticleSelection::InitControlObjects();
 135 }
 136
 137 THnSparse* AliDxHFEParticleSelectionD0::DefineTHnSparse()
 138 {
 139   //
 140   // Defines the THnSparse.
 141
 142   // here is the only place to change the dimension
 143   const int thnSize2 = 5;
 144   InitTHnSparseArray(thnSize2);
 145
 146   const double Pi=TMath::Pi();
 147   TString name;
 148   name.Form("%s info", GetName());
 149
 150   //                                 0     1     2       3         4
 151   //                                 Pt   Phi   Ptbin  D0InvMass  Eta
 152   int         thnBins [thnSize2] = {1000, 200,  15,     200,     500 };
 153   double      thnMin  [thnSize2] = {  0,    0,   0,    1.5648,   -1. };
 154   double      thnMax  [thnSize2] = { 100, 2*Pi, 14,    2.1648,    1. };
 155   const char* thnNames[thnSize2] = {"Pt", "Phi","Ptbin","D0InvMass","Eta"};
 156
 157   return CreateControlTHnSparse(name,thnSize2,thnBins,thnMin,thnMax,thnNames);
 158 }
 159
 160 int AliDxHFEParticleSelectionD0::FillParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
 161 {
 162   // fill the data array from the particle data
 163   if (!data) return -EINVAL;
 164   AliAODRecoDecayHF2Prong* track=dynamic_cast<AliAODRecoDecayHF2Prong*>(p);
 165   if (!track) return -ENODATA;
 166   int i=0;
 167   if (dimension!=GetDimTHnSparse()) {
 168     // TODO: think about filling only the available data and throwing a warning
 169     return -ENOSPC;
 170   }
 171   data[i++]=track->Pt();
 172   data[i++]=track->Phi();
 173   data[i++]=fPtBin;
 174   data[i++]=fD0InvMass;
 175   data[i++]=track->Eta();
 176
 177   return i;
 178 }
 179
 180 int AliDxHFEParticleSelectionD0::InitControlObjectsDaughters(TString name, int daughter)
 181 {
 182   // Setting up Control objects for the daughters.
 183   // Move to ParticleSelection??
 184   AliInfo("Setting up daughter THnSparse");
 185
 186   const int thnSize2 = 5;
 187   const double Pi=TMath::Pi();
 188   //                           0    1      2      3          4
 189   //                           Pt   Phi   Ptbin  D0InvMass  Eta
 190   int    thnBins[thnSize2] = { 1000,  200, 21,     200,     500};
 191   double thnMin [thnSize2] = {    0,    0,  0,    1.5648,   -1.};
 192   double thnMax [thnSize2] = {  100, 2*Pi, 20,    2.1648,    1.};
 193
 194   std::auto_ptr<THnSparseF> DaughterProperties(new THnSparseF(name, name, thnSize2, thnBins, thnMin, thnMax));
 195
 196   if (DaughterProperties.get()==NULL) {
 197     return -ENOMEM;
 198   }
 199
 200   for(int iLabel=0; iLabel< 5;iLabel++)
 201     DaughterProperties->GetAxis(iLabel)->SetTitle(fgkDgTrackControlBinNames[iLabel]);
 202
 203   if(daughter==0){
 204     fD0Daughter0=DaughterProperties.release();
 205     AddControlObject(fD0Daughter0);
 206   }
 207
 208   if(daughter==1){
 209     fD0Daughter1=DaughterProperties.release();
 210     AddControlObject(fD0Daughter1);
 211   }
 212   return 0;
 213 }
 214
 215 int AliDxHFEParticleSelectionD0::HistogramParticleProperties(AliVParticle* p, int selectionCode)
 216 {
 217
 218   /// histogram particle properties
 219   if (!p) return -EINVAL;
 220
 221   // fill the common histograms
 222   AliDxHFEParticleSelection::HistogramParticleProperties(p, selectionCode);
 223
 224   // no daughters to fill if 0 (= no candidate)
 225   if (selectionCode==0) return 0;
 226
 227   AliAODRecoDecayHF2Prong* part=dynamic_cast<AliAODRecoDecayHF2Prong*>(p);
 228
 229   if(!part) return 0;
 230   // Convention: 1. daughter is postive track, 2. = negative
 231   AliAODTrack *prongpos=(AliAODTrack*)part->GetDaughter(0);
 232   AliAODTrack *prongneg=(AliAODTrack*)part->GetDaughter(1);
 233
 234   if(!prongpos || !prongneg) {
 235     return 0;
 236   }
 237
 238   fD0InvMass= part->InvMassD0();
 239   fPtBin=fCuts->PtBin(part->Pt());
 240
 241   // TODO: avoid repeated allocation of the arrays
 242   Double_t KProperties[]={prongneg->Pt(),prongneg->Phi(),(Double_t)fPtBin, fD0InvMass,prongneg->Eta()};
 243   Double_t piProperties[]={prongpos->Pt(),prongpos->Phi(),(Double_t)fPtBin,fD0InvMass,prongpos->Eta()};
 244
 245
 246   // Fills only for D0 or both..
 247   if ((selectionCode==1 || selectionCode==3) && fFillOnlyD0D0bar<2) {
 248
 249     if(fD0Properties && ParticleProperties()) {
 250       memset(ParticleProperties(), 0, GetDimTHnSparse()*sizeof(ParticleProperties()[0]));
 251       FillParticleProperties(p, ParticleProperties(), GetDimTHnSparse());
 252       fD0Properties->Fill(ParticleProperties());
 253     }
 254     if(fD0Daughter0) fD0Daughter0->Fill(piProperties);
 255     if(fD0Daughter1) fD0Daughter1->Fill(KProperties);
 256   }
 257   // Checks for D0bar (or hypothesis both)
 258   if ((selectionCode==2 || selectionCode==3) && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) {
 259     // Set the fD0InvMass to InvMassD0bar instead..
 260     fD0InvMass= part->InvMassD0bar();
 261     if(fD0Properties && ParticleProperties()) {
 262       memset(ParticleProperties(), 0, GetDimTHnSparse()*sizeof(ParticleProperties()[0]));
 263       FillParticleProperties(p, ParticleProperties(), GetDimTHnSparse());
 264       fD0Properties->Fill(ParticleProperties());
 265     }
 266     if(fD0Daughter0) fD0Daughter0->Fill(piProperties);
 267     if(fD0Daughter1) fD0Daughter1->Fill(KProperties);
 268     //reset value to InvMassD0 for when CreateParticle() is called
 269     fD0InvMass= part->InvMassD0();
 270
 271   }
 272   return 0;
 273 }
 274
 275 TObjArray* AliDxHFEParticleSelectionD0::Select(TObjArray* pTracks, const AliVEvent *pEvent)
 276 {
 277   /// create selection, array contains only pointers but does not own the objects
 278   /// object array needs to be deleted by caller
 279   if (!pTracks) return NULL;
 280   TObjArray* selectedTracks=new TObjArray;
 281   if (!selectedTracks) return NULL;
 282   selectedTracks->SetOwner(kFALSE);
 283   TIter itrack(pTracks);
 284   TObject* pObj=NULL;
 285   while ((pObj=itrack())!=NULL) {
 286     AliVParticle* track=dynamic_cast<AliVParticle*>(pObj);
 287     if (!track) continue;
 288     int selectionCode=IsSelected(track,pEvent);
 289     HistogramParticleProperties(track, selectionCode);
 290
 291     // Add track if it is either defined as D0(selectionCode==1) or both
 292     // D0bar and a D0 (selectionCode==3)
 293     if ((selectionCode==1 || selectionCode==3) && fFillOnlyD0D0bar<2)
 294       selectedTracks->Add(CreateParticle(track));
 295
 296     // Add track if it is either defined as D0bar(selectionCode==2) or both
 297     // D0bar and a D0 (selectionCode==3)
 298     if ((selectionCode==2 || selectionCode==3) && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)){
 299       AliAODRecoDecayHF2Prong* prong=dynamic_cast<AliAODRecoDecayHF2Prong*>(track);
 300       fD0InvMass=prong?prong->InvMassD0bar():0.;
 301       selectedTracks->Add(CreateParticle(track));
 302     }
 303   }
 304
 305   HistogramEventProperties(AliDxHFEParticleSelection::kHistoNrTracksPrEvent,selectedTracks->GetEntries());
 306   return selectedTracks;
 307 }
 308
 309 int AliDxHFEParticleSelectionD0::IsSelected(AliVParticle* p, const AliVEvent* pEvent)
 310 {
 311   /// TODO: implement specific selection of D0 candidates
 312   /// Could also return values based on where where selection "failed
 313   /// Selected. Return 0 (none), 1(D0), 2(D0bar) or 3 (both)
 314
 315   int selectionCode=0;
 316
 317   AliAODRecoDecayHF2Prong *d0 = dynamic_cast<AliAODRecoDecayHF2Prong*>(p);
 318   if (!d0) return 0;
 319   if(d0->GetSelectionMap()) if(!d0->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
 320       AliDebug(1,"Skip D0 from Dstar");
 321       ((TH1D*)fHistoList->FindObject("fWhichCutD0"))->Fill(kDstar);
 322
 323       return 0; //skip the D0 from Dstar
 324     }
 325
 326   // TODO: the cuts instance should be const but the function definition of
 327   // AliRDHFCuts::IsSelected does not allow this
 328   AliRDHFCuts* cuts=const_cast<AliRDHFCuts*>(fCuts);
 329   if (!cuts) {
 330     selectionCode=0;
 331   }
 332   else if(cuts->IsInFiducialAcceptance(d0->Pt(),d0->Y(421)) ) {
 333
 334     // TODO: the aod pointer should also be const but the function definition of
 335     // AliRDHFCuts::IsSelected does not allow this
 336     AliAODEvent* aod=NULL;
 337     if (pEvent) aod=dynamic_cast<AliAODEvent*>(const_cast<AliVEvent*>(pEvent));
 338
 339     //TODO: Should add fSystem for PbPb    if(fSys==0){
 340     if(cuts->IsSelected(d0,AliRDHFCuts::kTracks,aod))       ((TH1D*)fHistoList->FindObject("fWhichCutD0"))->Fill(kCandSelTrack);
 341
 342     Int_t ptbin=cuts->PtBin(d0->Pt());
 343     if(ptbin==-1) {
 344       ((TH1D*)fHistoList->FindObject("fWhichCutD0"))->Fill(kNegPtbin);
 345       AliDebug(1,"Pt out of bounds");
 346       return 0;
 347     } //out of bounds
 348
 349     // Selected. Return 0 (none), 1 (D0), 2 (D0bar) or 3 (both)
 350     selectionCode=cuts->IsSelected(d0,AliRDHFCuts::kAll,aod);
 351     if(selectionCode==0)
 352       ((TH1D*)fHistoList->FindObject("fWhichCutD0"))->Fill(kSelected0);
 353
 354     if(selectionCode==1)
 355       ((TH1D*)fHistoList->FindObject("fWhichCutD0"))->Fill(kSelectedD0);
 356
 357     if(selectionCode==2)
 358       ((TH1D*)fHistoList->FindObject("fWhichCutD0"))->Fill(kSelectedD0bar);
 359
 360     if(selectionCode==3)
 361       ((TH1D*)fHistoList->FindObject("fWhichCutD0"))->Fill(kSelectedboth);
 362
 363     AliDebug(1,Form("Candidate is %d \n", selectionCode));
 364
 365     // check daughters before calling as there is unchecked code in
 366     // AliAODRecoDecayHF::HasBadDaughters called
 367     TObject* o=NULL;
 368     if (!((o=d0->GetDaughter(0))!=NULL && dynamic_cast<AliAODTrack*>(o)!=NULL &&
 369           (o=d0->GetDaughter(1))!=NULL && dynamic_cast<AliAODTrack*>(o)!=NULL)) {
 370       ((TH1D*)fHistoList->FindObject("fWhichCutD0"))->Fill(kNoDaugthers);
 371       AliDebug(1,"at least one daughter not found!");
 372
 373     }
 374   }
 375   else{
 376     ((TH1D*)fHistoList->FindObject("fWhichCutD0"))->Fill(kIsInFinducialAcceptance);
 377   }
 378
 379   return selectionCode;
 380 }
 381
 382 void AliDxHFEParticleSelectionD0::SetCuts(TObject* cuts, int level)
 383 {
 384   /// set cuts objects
 385   if (level==kCutD0) {
 386     fCuts=dynamic_cast<AliRDHFCuts*>(cuts);
 387     if (!fCuts && cuts) {
 388       AliError(Form("cuts object is not of required type AliRDHFCuts but %s", cuts->ClassName()));
 389     }
 390     return;
 391   }
 392   if (level==kCutList){
 393     TList* CutList=dynamic_cast<TList*>(cuts);
 394     if (!CutList && cuts) {
 395       AliError(Form("cuts object is not of required type TList but %s", cuts->ClassName()));
 396     }
 397     else{
 398       TObject *obj=NULL;
 399       int iii=0;
 400       TIter next(CutList);
 401       while((obj = next())){
 402         iii++;
 403         if(iii==1) {
 404           fCuts=dynamic_cast<AliRDHFCuts*>(obj);
 405           if (!fCuts)
 406             AliError(Form("Cut object is not of required type AliRDHFCuts but %s", obj->ClassName()));
 407         }
 408       }
 409     }
 410     return;
 411   }
 412   return;
 413 }
 414
 415 int AliDxHFEParticleSelectionD0::ParseArguments(const char* arguments)
 416 {
 417   // parse arguments and set internal flags
 418   TString strArguments(arguments);
 419   auto_ptr<TObjArray> tokens(strArguments.Tokenize(" "));
 420   if (!tokens.get()) return 0;
 421
 422   AliInfo(strArguments);
 423   TIter next(tokens.get());
 424   TObject* token;
 425   while ((token=next())) {
 426     TString argument=token->GetName();
 427     if (argument.BeginsWith("fillD0scheme=")){
 428       argument.ReplaceAll("fillD0scheme=", "");
 429       if (argument.CompareTo("both")==0){ fFillOnlyD0D0bar=0;}
 430       else if (argument.CompareTo("D0")==0){ fFillOnlyD0D0bar=1;}
 431       else if (argument.CompareTo("D0bar")==0){ fFillOnlyD0D0bar=2;}
 432       else {
 433         AliWarning(Form("can not set D0 filling scheme, unknown parameter '%s'", argument.Data()));
 434         fFillOnlyD0D0bar=0;
 435       }
 436       continue;
 437     }
 438     // forwarding of single argument works, unless key-option pairs separated
 439     // by blanks are introduced
 440     AliDxHFEParticleSelection::ParseArguments(argument);
 441   }
 442
 443   return 0;
 444 }
 445
 446 AliVParticle *AliDxHFEParticleSelectionD0::CreateParticle(AliVParticle* track)
 447 {
 448   //
 449   // Creates object containing only the variables needed for correlation
 450   //
 451
 452   AliReducedParticle *part = new AliReducedParticle(track->Eta(), track->Phi(), track->Pt(),fD0InvMass,fPtBin);
 453
 454   return part;
 455
 456 }