PWG3/vertexingHF/AliAnalysisTaskSED0Mass.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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 /////////////////////////////////////////////////////////////
  17 //
  18 // AliAnalysisTaskSE for D0 candidates invariant mass histogram
  19 // and comparison with the MC truth.
  20 //
  21 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
  22 // and Chiara Bianchin, chiara.bianchin@pd.infn.it
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <Riostream.h>
  26 #include <TClonesArray.h>
  27 #include <TNtuple.h>
  28 #include <TList.h>
  29 #include <TH1F.h>
  30
  31
  32 #include "AliAODEvent.h"
  33 #include "AliAODVertex.h"
  34 #include "AliAODTrack.h"
  35 #include "AliAODMCHeader.h"
  36 #include "AliAODMCParticle.h"
  37 #include "AliAODRecoDecayHF2Prong.h"
  38 #include "AliAODRecoCascadeHF.h"
  39 #include "AliAnalysisVertexingHF.h"
  40 #include "AliAnalysisTaskSE.h"
  41 #include "AliAnalysisTaskSED0Mass.h"
  42
  43 ClassImp(AliAnalysisTaskSED0Mass)
  44
  45
  46 //________________________________________________________________________
  47 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
  48 AliAnalysisTaskSE(),
  49 fOutputtight(0),
  50 fOutputloose(0),
  51 //fNtupleD0Cmp(0),
  52 /*
  53 fhistMass(0),
  54 fhistSgn(0),
  55 fhistBkg(0),
  56 */
  57 fVHFtight(0),
  58 fVHFloose(0),
  59 fNentries(0)
  60 {
  61   // Default constructor
  62 }
  63
  64 //________________________________________________________________________
  65 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name):
  66 AliAnalysisTaskSE(name),
  67 fOutputtight(0),
  68 fOutputloose(0),
  69 //fNtupleD0Cmp(0),
  70 /*
  71 fhistMass(0),
  72 fhistSgn(0),
  73 fhistBkg(0),
  74 */
  75 fVHFtight(0),
  76 fVHFloose(0),
  77 fNentries(0)
  78 {
  79   // Default constructor
  80
  81   // Output slot #1 writes into a TList container
  82   DefineOutput(1,TList::Class());  //My private output
  83   // Output slot #2 writes into a TList container
  84   DefineOutput(2,TList::Class());  //My private output
  85   // Output slot #3 writes into a TH1F container
  86   DefineOutput(3,TH1F::Class());  //My private output
  87 }
  88
  89 //________________________________________________________________________
  90 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
  91 {
  92   // Destructor
  93   /*
  94   if (fhistMass){
  95     delete fhistMass;
  96     fhistMass=0;
  97   }
  98
  99   if (fhistSgn){
 100     delete fhistSgn;
 101     fhistSgn=0;
 102   }
 103
 104   if (fhistBkg){
 105     delete fhistBkg;
 106     fhistBkg=0;
 107   }
 108   */
 109   if (fOutputtight) {
 110     delete fOutputtight;
 111     fOutputtight = 0;
 112   }
 113   if (fVHFtight) {
 114     delete fVHFtight;
 115     fVHFtight = 0;
 116   }
 117   if (fOutputloose) {
 118     delete fOutputloose;
 119     fOutputloose = 0;
 120   }
 121   if (fVHFloose) {
 122     delete fVHFloose;
 123     fVHFloose = 0;
 124   }
 125
 126 }
 127
 128 //________________________________________________________________________
 129 void AliAnalysisTaskSED0Mass::Init()
 130 {
 131   // Initialization
 132
 133   if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
 134
 135   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/ConfigVertexingHF.C");
 136
 137   // 2 sets of dedidcated cuts -- defined in UserExec
 138   fVHFtight = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
 139   fVHFloose = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
 140
 141   return;
 142 }
 143
 144 //________________________________________________________________________
 145 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
 146 {
 147
 148   // Create the output container
 149   //
 150   if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
 151
 152   // Several histograms are more conveniently managed in a TList
 153   fOutputtight = new TList();
 154   fOutputtight->SetOwner();
 155   fOutputtight->SetName("listtight");
 156
 157   fOutputloose = new TList();
 158   fOutputloose->SetOwner();
 159   fOutputloose->SetName("listloose");
 160
 161   const Int_t nhist=4;
 162
 163   TString nameMass=" ", nameSgn=" ", nameBkg=" ", nameRfl=" ";
 164
 165   for(Int_t i=0;i<nhist;i++){
 166     nameMass="histMass_";
 167     nameMass+=i+1;
 168     nameSgn="histSgn_";
 169     nameSgn+=i+1;
 170     nameBkg="histBkg_";
 171     nameBkg+=i+1;
 172     nameRfl="histRfl_";
 173     nameRfl+=i+1;
 174
 175     TH1F* tmpM = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
 176     tmpM->Sumw2();
 177
 178     TH1F* tmpS = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
 179     tmpS->Sumw2();
 180
 181     TH1F* tmpB = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
 182     tmpB->Sumw2();
 183
 184     //Reflection: histo filled with D0 which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 185     TH1F* tmpR = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
 186     tmpR->Sumw2();
 187
 188   /*
 189     tmpM->SetName(nameMass);
 190     tmpS->SetName(nameSgn);
 191     tmpB->SetName(nameBkg);
 192     tmpR->SetName(nameRfl);
 193   */
 194   //  printf("Created histograms %s\t%s\t%s\t%s\n",tmpM->GetName(),tmpS->GetName(),tmpB->GetName(),tmpR->GetName());
 195
 196     fOutputtight->Add(tmpM);
 197     fOutputtight->Add(tmpS);
 198     fOutputtight->Add(tmpB);
 199     fOutputtight->Add(tmpR);
 200
 201     fOutputloose->Add(tmpM);
 202     fOutputloose->Add(tmpS);
 203     fOutputloose->Add(tmpB);
 204     fOutputloose->Add(tmpR);
 205
 206 //   delete tmpM;
 207 //   delete tmpS;
 208 //   delete tmpB;
 209 //   delete tmpR;
 210
 211   }
 212
 213
 214   //fOutputloose->ls();
 215   //fOutputtight->ls();
 216
 217   fNentries=new TH1F("nentries", "Look at the number of entries! = number of AODs", 2,1.,2.);
 218
 219   return;
 220 }
 221
 222 //________________________________________________________________________
 223 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
 224 {
 225   // Execute analysis for current event:
 226   // heavy flavor candidates association to MC truth
 227   //cout<<"I'm in UserExec"<<endl;
 228   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 229
 230   // load D0->Kpi candidates
 231   TClonesArray *inputArrayD0toKpi =
 232     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 233   if(!inputArrayD0toKpi) {
 234     printf("AliAnalysisTaskSECompareHFpt::UserExec: D0toKpi branch not found!\n");
 235     return;
 236   }
 237
 238   // AOD primary vertex
 239   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 240   //vtx1->Print();
 241
 242   // load MC particles
 243   TClonesArray *mcArray =
 244     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 245   if(!mcArray) {
 246     printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
 247     return;
 248   }
 249
 250   // load MC header
 251   AliAODMCHeader *mcHeader =
 252     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 253   if(!mcHeader) {
 254     printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
 255     return;
 256   }
 257
 258
 259   //histogram filled with 1 for every AOD
 260   fNentries->Fill(1);
 261   PostData(3,fNentries);
 262   // loop over D0->Kpi candidates
 263   Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast();
 264   printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 265
 266   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 267
 268     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayD0toKpi->UncheckedAt(iD0toKpi);
 269     Bool_t unsetvtx=kFALSE;
 270     if(!d->GetOwnPrimaryVtx()) {
 271       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 272       unsetvtx=kTRUE;
 273     }
 274
 275     //  cout<<"Cuts applied: ";
 276     //     Double_t *cutsapplied=(Double_t*)fVHF->GetD0toKpiCuts();
 277     //     for (Int_t i=0;i<9;i++){
 278     //       printf(cutsapplied[i]\t);
 279     //     }
 280     //    cout<<endl;
 281
 282     //cuts order
 283 //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
 284 //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
 285 //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
 286 //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
 287 //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
 288 //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
 289 //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
 290 //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
 291 //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
 292
 293
 294     Double_t pt = d->Pt();
 295     Int_t ptbin=0;
 296
 297     //cout<<"P_t = "<<pt<<endl;
 298     if (pt>0. && pt<=1.) {
 299       ptbin=1;
 300       fVHFtight->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.0003,0.7);
 301       fVHFloose->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.00025,0.7);
 302 //       fVHFtight->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.0002,0.7);
 303 //       fVHFloose->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,1,1,-0.00015,0.5);
 304       //printf("I'm in the bin %d\n",ptbin);
 305
 306     }
 307     else {
 308       if(pt>1. && pt<=3.) {
 309         ptbin=2;
 310         fVHFtight->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.0003,0.9);
 311         fVHFloose->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,1,1,-0.00025,0.8);
 312         //printf("I'm in the bin %d\n",ptbin);
 313
 314       }
 315        if(pt>3. && pt<=5.){
 316         ptbin=3;
 317         fVHFtight->SetD0toKpiCuts(0.7,0.015,0.8,0.7,0.7,0.05,0.05,-0.0002,0.9);
 318         fVHFloose->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00015,0.8);
 319         //printf("I'm in the bin %d\n",ptbin);
 320
 321       } else{
 322         ptbin=4;
 323         fVHFtight->SetD0toKpiCuts(0.7,0.015,0.8,0.7,0.7,0.05,0.05,-0.0002,0.95);
 324         fVHFloose->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00015,0.9);
 325       }//if(pt>5)
 326     }
 327     //printf("I'm in the bin %d\n",ptbin);
 328     //old
 329     //fVHF->SetD0toKpiCuts(0.7,0.03,0.8,0.06,0.06,0.05,0.05,-0.0002,0.6); //2.p-p vertex reconstructed
 330
 331
 332     FillHists(ptbin,d,mcArray,fVHFtight,fOutputtight);
 333     FillHists(ptbin,d,mcArray,fVHFloose,fOutputloose);
 334     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 335
 336   }
 337
 338
 339
 340   // Post the data
 341   PostData(1,fOutputtight);
 342   PostData(2,fOutputloose);
 343
 344   return;
 345 }
 346 //____________________________________________________________________________*
 347 void AliAnalysisTaskSED0Mass::FillHists(Int_t ptbin, AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAnalysisVertexingHF *vhf, TList *listout){
 348   //
 349   // function used in UserExec:
 350   //
 351   Int_t okD0=0,okD0bar=0;
 352
 353   if(part->SelectD0(vhf->GetD0toKpiCuts(),okD0,okD0bar)) {//selected
 354     Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
 355     //printf("SELECTED\n");
 356     Int_t labD0 = part->MatchToMC(421,arrMC); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
 357     //printf("labD0 %d",labD0);
 358
 359     TString fillthis="";
 360     TString fillrfl= "";
 361
 362     if (okD0==1) {
 363       fillthis="histMass_";
 364       fillthis+=ptbin;
 365       //cout<<"Filling "<<fillthis<<endl;
 366
 367       //printf("Fill mass with D0");
 368       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 369       if(labD0>=0) {
 370         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 371         Int_t pdgD0 = partD0->GetPdgCode();
 372         //cout<<"pdg = "<<pdgD0<<endl;
 373         fillthis="histSgn_";
 374         fillrfl="histRfl_";
 375         fillthis+=ptbin;
 376         fillrfl+=ptbin;
 377         if (pdgD0==421){ //D0
 378           //cout<<"Fill S with D0"<<endl;
 379           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 380
 381         } else{ //it was a D0bar
 382           ((TH1F*)(listout->FindObject(fillrfl)))->Fill(invmassD0);
 383         }
 384       } else {//background
 385         fillthis="histBkg_";
 386         fillthis+=ptbin;
 387         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 388       }
 389     }
 390     if (okD0bar==1) {
 391       fillthis="histMass_";
 392       fillthis+=ptbin;
 393       //printf("Fill mass with D0bar");
 394       ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
 395
 396       if(labD0>=0) {
 397         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 398         Int_t pdgD0 = partD0->GetPdgCode();
 399         //cout<<" pdg = "<<pdgD0<<endl;
 400         if (pdgD0==-421){ //D0bar
 401           fillthis="histSgn_";
 402           fillthis+=ptbin;
 403           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 404
 405         } else{
 406           ((TH1F*)(listout->FindObject(fillrfl)))->Fill(invmassD0bar);
 407
 408         }
 409       } else {//background
 410         fillthis="histBkg_";
 411         fillthis+=ptbin;
 412         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 413       }
 414     }
 415
 416
 417   } //else cout<<"NOT SELECTED"<<endl;
 418
 419 }
 420 //________________________________________________________________________
 421 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
 422 {
 423   // Terminate analysis
 424   //
 425   if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
 426
 427   fOutputtight = dynamic_cast<TList*> (GetOutputData(1));
 428   if (!fOutputtight) {
 429     printf("ERROR: fOutput not available\n");
 430     return;
 431   }
 432   fOutputloose = dynamic_cast<TList*> (GetOutputData(2));
 433   if (!fOutputloose) {
 434     printf("ERROR: fOutput not available\n");
 435     return;
 436   }
 437
 438
 439   return;
 440 }
 441