PWGHF/vertexingHF/AliAnalysisTaskSESignificance.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, 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 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // AliAnalysisTaskSESignificane to calculate effects on
  21 // significance of D mesons  cut
  22 // Authors: G. Ortona, ortona@to.infn.it
  23 // F. Prino, prino@to.infn.it
  24 // Renu Bala, bala@to.infn.it
  25 // Chiara Bianchin, cbianchi@pd.infn.it
  26 /////////////////////////////////////////////////////////////
  27
  28 #include <Riostream.h>
  29 #include <TClonesArray.h>
  30 #include <TCanvas.h>
  31 #include <TList.h>
  32 #include <TFile.h>
  33 #include <TString.h>
  34 #include <TH1F.h>
  35 #include <TDatabasePDG.h>
  36
  37 #include <AliLog.h>
  38 #include "AliAnalysisManager.h"
  39 #include "AliAODHandler.h"
  40 #include "AliAODEvent.h"
  41 #include "AliAODVertex.h"
  42 #include "AliAODTrack.h"
  43 #include "AliAODMCHeader.h"
  44 #include "AliAODMCParticle.h"
  45 #include "AliAODRecoDecayHF3Prong.h"
  46 #include "AliAODRecoDecayHF.h"
  47 #include "AliAODRecoDecayHF2Prong.h"
  48 #include "AliAODRecoDecayHF4Prong.h"
  49 #include "AliAODRecoCascadeHF.h"
  50
  51 #include "AliAnalysisTaskSE.h"
  52 #include "AliRDHFCutsDplustoKpipi.h"
  53 #include "AliRDHFCutsD0toKpipipi.h"
  54 #include "AliRDHFCutsDstoKKpi.h"
  55 #include "AliRDHFCutsDStartoKpipi.h"
  56 #include "AliRDHFCutsD0toKpi.h"
  57 #include "AliRDHFCutsLctopKpi.h"
  58 #include "AliMultiDimVector.h"
  59
  60 #include "AliAnalysisTaskSESignificance.h"
  61
  62 using std::cout;
  63 using std::endl;
  64
  65 ClassImp(AliAnalysisTaskSESignificance)
  66
  67
  68 //________________________________________________________________________
  69 AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance():
  70   AliAnalysisTaskSE(),
  71   fOutput(0),
  72   fCutList(0),
  73   fHistNEvents(0),
  74   fUpmasslimit(1.965),
  75   fLowmasslimit(1.765),
  76   fRDCuts(0),
  77   fNPtBins(0),
  78   fReadMC(kFALSE),
  79   fUseSelBit(kFALSE),
  80   fBFeedDown(kBoth),
  81   fDecChannel(0),
  82   fPDGmother(0),
  83   fNProngs(0),
  84   fBranchName(""),
  85   fSelectionlevel(0),
  86   fNVars(0),
  87   fNBins(100),
  88   fPartOrAndAntiPart(0),
  89   fDsChannel(0)
  90 {
  91   // Default constructor
  92   SetPDGCodes();
  93   SetDsChannel(kPhi);
  94
  95   for(Int_t i=0;i<4;i++) fPDGdaughters[i]=0;
  96   for(Int_t i=0;i<2;i++) {fPDGD0ToKpi[i]=0;fPDGDStarToD0pi[i]=0;}
  97   for(Int_t i=0;i<kMaxCutVar;i++) fVars[i]=0.;
  98   for(Int_t i=0;i<kMaxNHist;i++) {
  99     fMassHist[i]=0;
 100     fSigHist[i]=0;
 101     fBkgHist[i]=0;
 102     fRflHist[i]=0;
 103   }
 104
 105 }
 106
 107 //________________________________________________________________________
 108 AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance(const char *name, TList* listMDV,AliRDHFCuts *rdCuts,Int_t decaychannel,Int_t selectionlevel):
 109   AliAnalysisTaskSE(name),
 110   fOutput(0),
 111   fCutList(listMDV),
 112   fHistNEvents(0),
 113   fUpmasslimit(0),
 114   fLowmasslimit(0),
 115   fRDCuts(rdCuts),
 116   fNPtBins(0),
 117   fReadMC(kFALSE),
 118   fUseSelBit(kFALSE),
 119   fBFeedDown(kBoth),
 120   fDecChannel(decaychannel),
 121   fPDGmother(0),
 122   fNProngs(0),
 123   fBranchName(""),
 124   fSelectionlevel(selectionlevel),
 125   fNVars(0),
 126   fNBins(100),
 127   fPartOrAndAntiPart(0),
 128   fDsChannel(0)
 129 {
 130
 131   SetPDGCodes();
 132   SetDsChannel(kPhi);
 133   if (fDecChannel!=2) SetMassLimits(0.15,fPDGmother); //check range
 134   else {
 135     Float_t min = 0.13;
 136     Float_t max = 0.19;
 137     SetMassLimits(min, max);
 138   }
 139   fNPtBins=fRDCuts->GetNPtBins();
 140
 141   fNVars=fRDCuts->GetNVarsForOpt();
 142   if(fNVars>kMaxCutVar) AliFatal(Form("Too large number of cut variables, maximum is %d",kMaxCutVar));
 143
 144   if(fDebug>1)fRDCuts->PrintAll();
 145    // Output slot #1 writes into a TList container
 146   DefineOutput(1,TList::Class());  //My private output
 147   DefineOutput(2,TList::Class());
 148   DefineOutput(3,AliRDHFCuts::Class()); //class of the cuts
 149   CheckConsistency();
 150 }
 151
 152  //________________________________________________________________________
 153 AliAnalysisTaskSESignificance::~AliAnalysisTaskSESignificance()
 154 {
 155   // Destructor
 156   if (fOutput) {
 157     delete fOutput;
 158     fOutput = 0;
 159   }
 160   if (fCutList) {
 161     delete fCutList;
 162     fCutList = 0;
 163   }
 164   if(fHistNEvents){
 165     delete fHistNEvents;
 166     fHistNEvents=0;
 167   }
 168 /*
 169   if(fRDCuts) {
 170     delete fRDCuts;
 171     fRDCuts= 0;
 172   }
 173 */
 174
 175 }
 176 //_________________________________________________________________
 177 void AliAnalysisTaskSESignificance::SetPDGCodes(){
 178   // sets channel dependent quantities
 179
 180   switch(fDecChannel){
 181   case 0:
 182     //D+
 183     fPDGmother=411;
 184     fNProngs=3;
 185     fPDGdaughters[0]=211;//pi
 186     fPDGdaughters[1]=321;//K
 187     fPDGdaughters[2]=211;//pi
 188     fPDGdaughters[3]=0; //empty
 189     fBranchName="Charm3Prong";
 190     break;
 191   case 1:
 192     //D0
 193     fPDGmother=421;
 194     fNProngs=2;
 195     fPDGdaughters[0]=211;//pi
 196     fPDGdaughters[1]=321;//K
 197     fPDGdaughters[2]=0; //empty
 198     fPDGdaughters[3]=0; //empty
 199     fBranchName="D0toKpi";
 200     break;
 201   case 2:
 202     //D*
 203     fPDGmother=413;
 204     fNProngs=3;
 205     fPDGdaughters[1]=211;//pi
 206     fPDGdaughters[0]=321;//K
 207     fPDGdaughters[2]=211;//pi (soft?)
 208     fPDGdaughters[3]=0; //empty
 209     fBranchName="Dstar";
 210     break;
 211   case 3:
 212     //Ds
 213     fPDGmother=431;
 214     fNProngs=3;
 215     fPDGdaughters[0]=321;//K
 216     fPDGdaughters[1]=321;//K
 217     fPDGdaughters[2]=211;//pi
 218     fPDGdaughters[3]=0; //empty
 219     fBranchName="Charm3Prong";
 220     break;
 221   case 4:
 222     //D0 in 4 prongs
 223     fPDGmother=421;
 224     fNProngs=4;
 225     fPDGdaughters[0]=321;
 226     fPDGdaughters[1]=211;
 227     fPDGdaughters[2]=211;
 228     fPDGdaughters[3]=211;
 229     fBranchName="Charm4Prong";
 230     break;
 231   case 5:
 232     //Lambda_c
 233     fPDGmother=4122;
 234     fNProngs=3;
 235     fPDGdaughters[0]=2212;//p
 236     fPDGdaughters[1]=321;//K
 237     fPDGdaughters[2]=211;//pi
 238     fPDGdaughters[3]=0; //empty
 239     fBranchName="Charm3Prong";
 240     break;
 241   }
 242 }
 243
 244 //_________________________________________________________________
 245 Bool_t AliAnalysisTaskSESignificance::CheckConsistency(){
 246
 247   Bool_t result = kTRUE;
 248
 249   const Int_t nvars=fRDCuts->GetNVars();//ForOpt();
 250   //Float_t *vars = new Float_t[nvars];
 251   Bool_t *varsforopt = fRDCuts->GetVarsForOpt();
 252   Bool_t *uppervars = fRDCuts->GetIsUpperCut();
 253   TString *names = fRDCuts->GetVarNames();
 254
 255   for(Int_t i=0;i<fNPtBins;i++){
 256     TString mdvname=Form("multiDimVectorPtBin%d",i);
 257     Int_t ic=0;
 258
 259     for(Int_t ivar=0;ivar<nvars;ivar++){
 260       if(varsforopt[ivar]){
 261         Float_t min = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMinLimit(ic);
 262         Float_t max = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMaxLimit(ic);
 263         if(min==max){
 264           AliFatal(Form("tight and loose cut for optimization variable number %d are the same in ptbin %d\n",ic,i));
 265           return kFALSE;
 266         }
 267         Bool_t lowermdv = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGreaterThan(ic);
 268         if(uppervars[ivar]&&lowermdv){
 269           AliFatal(Form("%s is declared as uppercut, but as been given tighter cut larger then loose cut in ptbin %d \n ---please check your cuts \n ",names[ivar].Data(),i));
 270           return kFALSE;
 271         }
 272         if(!uppervars[ivar]&&!lowermdv){
 273           AliFatal(Form("%s is declared as lower cut, but as been given tighter cut smaller then loose cut in ptbin %d \n ---please check your cuts \n",names[ivar].Data(),i));
 274           return kFALSE;
 275         }
 276         ic++;
 277       }
 278     }
 279   }
 280   return result;
 281 }
 282 //_________________________________________________________________
 283 void AliAnalysisTaskSESignificance::SetBFeedDown(FeedDownEnum flagB){
 284   if(fReadMC)fBFeedDown=flagB;
 285   else {
 286     AliInfo("B feed down not allowed without MC info\n");
 287     fBFeedDown=kBoth;
 288   }
 289 }
 290 //_________________________________________________________________
 291 void  AliAnalysisTaskSESignificance::SetMassLimits(Float_t range, Int_t pdg){
 292   Float_t mass=0;
 293   Int_t abspdg=TMath::Abs(pdg);
 294   mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
 295   fUpmasslimit = mass+range;
 296   fLowmasslimit = mass-range;
 297 }
 298 //_________________________________________________________________
 299 void  AliAnalysisTaskSESignificance::SetMassLimits(Float_t lowlimit, Float_t uplimit){
 300   if(uplimit>lowlimit)
 301     {
 302       fUpmasslimit = uplimit;
 303       fLowmasslimit = lowlimit;
 304     }
 305 }
 306
 307
 308
 309 //________________________________________________________________________
 310 void AliAnalysisTaskSESignificance::LocalInit()
 311 {
 312   // Initialization
 313
 314   if(fDebug > 1) printf("AnalysisTaskSESignificance::Init() \n");
 315
 316   switch(fDecChannel){
 317   case 0:
 318     {
 319       AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCuts)));
 320       // Post the data
 321       PostData(3,copycut);
 322     }
 323     break;
 324   case 1:
 325     {
 326       AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCuts)));
 327       // Post the data
 328       PostData(3,copycut);
 329     }
 330     break;
 331   case 2:
 332     {
 333       AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCuts)));
 334       // Post the data
 335       PostData(3,copycut);
 336     }
 337     break;
 338   case 3:
 339     {
 340       AliRDHFCutsDstoKKpi* copycut=new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fRDCuts)));
 341       // Post the data
 342       PostData(3,copycut);
 343     }
 344     break;
 345   case 4:
 346     {
 347       AliRDHFCutsD0toKpipipi* copycut=new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fRDCuts)));
 348       // Post the data
 349       PostData(3,copycut);
 350     }
 351     break;
 352   case 5:
 353     {
 354       AliRDHFCutsLctopKpi* copycut=new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fRDCuts)));
 355       // Post the data
 356       PostData(3,copycut);
 357     }
 358     break;
 359
 360   default:
 361     return;
 362   }
 363
 364   TList *mdvList =  new TList();
 365   mdvList->SetOwner();
 366   mdvList = fCutList;
 367
 368   PostData(2,mdvList);
 369
 370
 371 }
 372 //________________________________________________________________________
 373 void AliAnalysisTaskSESignificance::UserCreateOutputObjects()
 374 {
 375   // Create the output container
 376
 377   if(fDebug > 1) printf("AnalysisTaskSESignificance::UserCreateOutputObjects() \n");
 378
 379   // Several histograms are more conveniently managed in a TList
 380   fOutput = new TList();
 381   fOutput->SetOwner();
 382   fOutput->SetName("OutputHistos");
 383
 384   //same number of steps in each multiDimVectorPtBin%d !
 385   Int_t nHist=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
 386   cout<<"ncells = "<<nHist<<" n ptbins = "<<fNPtBins<<endl;
 387   nHist=nHist*fNPtBins;
 388   cout<<"Total = "<<nHist<<endl;
 389   for(Int_t i=0;i<nHist;i++){
 390
 391     TString hisname;
 392     TString signame;
 393     TString bkgname;
 394     TString rflname;
 395     TString title;
 396
 397     hisname.Form("hMass_%d",i);
 398     signame.Form("hSig_%d",i);
 399     bkgname.Form("hBkg_%d",i);
 400     rflname.Form("hRfl_%d",i);
 401
 402     title.Form("Invariant mass;M[GeV/c^{2}];Entries");
 403
 404     fMassHist[i]=new TH1F(hisname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
 405     fMassHist[i]->Sumw2();
 406     fOutput->Add(fMassHist[i]);
 407
 408     if(fReadMC){
 409       fSigHist[i]=new TH1F(signame.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
 410       fSigHist[i]->Sumw2();
 411       fOutput->Add(fSigHist[i]);
 412
 413       fBkgHist[i]=new TH1F(bkgname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
 414       fBkgHist[i]->Sumw2();
 415       fOutput->Add(fBkgHist[i]);
 416
 417       if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi){
 418         fRflHist[i]=new TH1F(rflname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
 419         fRflHist[i]->Sumw2();
 420         fOutput->Add(fRflHist[i]);
 421       }
 422     }
 423   }
 424
 425   fHistNEvents=new TH1F("fHistNEvents","Number of AODs scanned",8,-0.5,7.5);
 426   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 427   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvSelected (vtx)");
 428   fHistNEvents->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
 429   fHistNEvents->GetXaxis()->SetBinLabel(4,"nTotEntries Mass hists");
 430   fHistNEvents->GetXaxis()->SetBinLabel(5,"Pile-up Rej");
 431   fHistNEvents->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
 432   if(fReadMC){
 433     fHistNEvents->GetXaxis()->SetBinLabel(7,"MC Cand from c");
 434     fHistNEvents->GetXaxis()->SetBinLabel(8,"MC Cand from b");
 435   } else{
 436     fHistNEvents->GetXaxis()->SetBinLabel(7,"N candidates");
 437   }
 438   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 439   fOutput->Add(fHistNEvents);
 440
 441
 442   PostData(1,fOutput);
 443   return;
 444 }
 445
 446 //________________________________________________________________________
 447 void AliAnalysisTaskSESignificance::UserExec(Option_t */*option*/)
 448 {
 449   // Execute analysis for current event:
 450   // heavy flavor candidates association to MC truth
 451
 452   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 453   if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
 454   // Post the data already here
 455   PostData(1,fOutput);
 456   TClonesArray *arrayProng =0;
 457
 458   if(!aod && AODEvent() && IsStandardAOD()) {
 459     // In case there is an AOD handler writing a standard AOD, use the AOD
 460     // event in memory rather than the input (ESD) event.
 461     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 462     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 463     // have to taken from the AOD event hold by the AliAODExtension
 464     AliAODHandler* aodHandler = (AliAODHandler*)
 465       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 466     if(aodHandler->GetExtensions()) {
 467
 468       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 469       AliAODEvent *aodFromExt = ext->GetAOD();
 470       arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
 471
 472    }
 473   } else if(aod) {
 474     arrayProng=(TClonesArray*)aod->GetList()->FindObject(fBranchName.Data());
 475   }
 476   if(!aod || !arrayProng) {
 477     AliError("AliAnalysisTaskSESignificance::UserExec:Branch not found!\n");
 478     return;
 479   }
 480
 481   // fix for temporary bug in ESDfilter
 482   // the AODs with null vertex pointer didn't pass the PhysSel
 483   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 484   TClonesArray *arrayMC=0;
 485   AliAODMCHeader *mcHeader=0;
 486
 487   // load MC particles
 488   if(fReadMC){
 489
 490     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 491     if(!arrayMC) {
 492       AliError("AliAnalysisTaskSESignificance::UserExec:MC particles branch not found!\n");
 493       return;
 494     }
 495
 496     // load MC header
 497     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 498     if(!mcHeader) {
 499       AliError("AliAnalysisTaskSESignificance::UserExec:MC header branch not found!\n");
 500       return;
 501     }
 502   }
 503
 504   fHistNEvents->Fill(0); // count event
 505
 506   AliAODRecoDecayHF *d=0;
 507
 508
 509
 510   Int_t nHistpermv=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
 511   Int_t nProng = arrayProng->GetEntriesFast();
 512   if(fDebug>1) printf("Number of D2H: %d\n",nProng);
 513
 514   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
 515   TString trigclass=aod->GetFiredTriggerClasses();
 516   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(5);
 517
 518   if(fRDCuts->IsEventSelected(aod)) {
 519     fHistNEvents->Fill(1);
 520   }else{
 521     if(fRDCuts->GetWhyRejection()==1) // rejected for pileup
 522       fHistNEvents->Fill(4);
 523     return;
 524   }
 525
 526   for (Int_t iProng = 0; iProng < nProng; iProng++) {
 527     fHistNEvents->Fill(6);
 528     d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
 529
 530     AliAODRecoCascadeHF* DStarToD0pi = NULL;
 531     AliAODRecoDecayHF2Prong* D0Particle = NULL;
 532     fPDGDStarToD0pi[0] = 421; fPDGDStarToD0pi[1] = 211;
 533     fPDGD0ToKpi[0] = 321; fPDGD0ToKpi[1] = 211;
 534
 535     Bool_t isSelBit=kTRUE;
 536     if(fUseSelBit){
 537       if(fDecChannel==0) {
 538         isSelBit=d->HasSelectionBit(AliRDHFCuts::kDplusCuts);
 539       }else{
 540         if(fDecChannel==1) {
 541           isSelBit=d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts);
 542         }else{
 543           if(fDecChannel==2) {
 544             isSelBit=d->HasSelectionBit(AliRDHFCuts::kDstarCuts);
 545           }else{
 546             if(fDecChannel==3) {
 547               isSelBit=d->HasSelectionBit(AliRDHFCuts::kDsCuts);
 548             }else{
 549               if(fDecChannel==5) isSelBit=d->HasSelectionBit(AliRDHFCuts::kLcCuts);
 550             }
 551           }
 552         }
 553       }
 554     }
 555     if(!isSelBit) continue;
 556
 557     if (fDecChannel==2) {
 558       DStarToD0pi = (AliAODRecoCascadeHF*)arrayProng->At(iProng);
 559       if (!DStarToD0pi->GetSecondaryVtx()) continue;
 560       D0Particle = (AliAODRecoDecayHF2Prong*)DStarToD0pi->Get2Prong();
 561       if (!D0Particle) continue;
 562     }
 563
 564     Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(fPDGmother));
 565     Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel,aod);
 566
 567     if(fReadMC && fBFeedDown!=kBoth && isSelected){
 568       Int_t labD = d->MatchToMC(fPDGmother,arrayMC,fNProngs,fPDGdaughters);
 569       if(labD>=0){
 570         AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
 571         Int_t pdgGranma = CheckOrigin(partD, arrayMC);
 572         Int_t abspdgGranma = TMath::Abs(pdgGranma);
 573         if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
 574           //feed down particle
 575           AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
 576           fHistNEvents->Fill(7);
 577           if(fBFeedDown==kCharmOnly) isSelected=kFALSE; //from beauty
 578         }
 579         else {
 580           //prompt particle
 581           fHistNEvents->Fill(6);
 582           if(fBFeedDown==kBeautyOnly)isSelected=kFALSE;
 583         }
 584       }
 585     }
 586
 587     if(isSelected&&isFidAcc) {
 588       fHistNEvents->Fill(2); // count selected with loosest cuts
 589       if(fDebug>1) printf("+++++++Is Selected\n");
 590
 591       Int_t nVals=0;
 592       if(fDecChannel==3) SetPDGdaughterDstoKKpi();
 593       fRDCuts->GetCutVarsForOpt(d,fVars,fNVars,fPDGdaughters,aod);
 594       Int_t ptbin=fRDCuts->PtBin(d->Pt());
 595       if(ptbin==-1) continue;
 596       TString mdvname=Form("multiDimVectorPtBin%d",ptbin);
 597       AliMultiDimVector* muvec=(AliMultiDimVector*)fCutList->FindObject(mdvname.Data());
 598
 599       ULong64_t *addresses = muvec->GetGlobalAddressesAboveCuts(fVars,(Float_t)d->Pt(),nVals);
 600       if(fDebug>1)printf("nvals = %d\n",nVals);
 601       for(Int_t ivals=0;ivals<nVals;ivals++){
 602         if(addresses[ivals]>=muvec->GetNTotCells()){
 603           if (fDebug>1) printf("Overflow!!\n");
 604           delete [] addresses;
 605           return;
 606         }
 607
 608         fHistNEvents->Fill(3);
 609
 610         //fill the histograms with the appropriate method
 611         switch (fDecChannel){
 612         case 0:
 613           FillDplus(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
 614           break;
 615         case 1:
 616           FillD02p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
 617           break;
 618         case 2:
 619           FillDstar(DStarToD0pi,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
 620           break;
 621         case 3:
 622           if(isSelected&1){
 623             FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected,1);
 624           }
 625           break;
 626         case 4:
 627           FillD04p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
 628           break;
 629         case 5:
 630           FillLambdac(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
 631           break;
 632         default:
 633           break;
 634         }
 635
 636       }
 637
 638       if (fDecChannel==3 && isSelected&2){
 639         SetPDGdaughterDstopiKK();
 640         nVals=0;
 641         fRDCuts->GetCutVarsForOpt(d,fVars,fNVars,fPDGdaughters,aod);
 642         delete [] addresses;
 643         addresses = muvec->GetGlobalAddressesAboveCuts(fVars,(Float_t)d->Pt(),nVals);
 644         if(fDebug>1)printf("nvals = %d\n",nVals);
 645         for(Int_t ivals=0;ivals<nVals;ivals++){
 646           if(addresses[ivals]>=muvec->GetNTotCells()){
 647             if (fDebug>1) printf("Overflow!!\n");
 648             delete [] addresses;
 649             return;
 650           }
 651           FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected,0);
 652
 653         }
 654
 655       }
 656
 657       delete [] addresses;
 658     }// end if selected
 659
 660   }
 661
 662
 663   PostData(1,fOutput);
 664   return;
 665 }
 666
 667 //***************************************************************************
 668
 669 // Methods used in the UserExec
 670
 671
 672 //********************************************************************************************
 673
 674 //Methods to fill the istograms with MC information, one for each candidate
 675 //NB: the implementation for each candidate is responsibility of the corresponding developer
 676
 677 void AliAnalysisTaskSESignificance::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
 678     //D+ channel
 679   if(!isSel){
 680     AliError("Candidate not selected\n");
 681     return;
 682   }
 683
 684   if(fPartOrAndAntiPart*d->GetCharge()<0)return;
 685   Int_t pdgdaughters[3] = {211,321,211};
 686   Double_t mass=d->InvMass(3,(UInt_t*)pdgdaughters);
 687
 688   fMassHist[index]->Fill(mass);
 689
 690
 691   if(fReadMC){
 692     Int_t lab=-1;
 693     lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
 694     if(lab>=0){ //signal
 695       fSigHist[index]->Fill(mass);
 696     } else{ //background
 697       fBkgHist[index]->Fill(mass);
 698     }
 699   }
 700 }
 701
 702
 703 void AliAnalysisTaskSESignificance::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
 704
 705   //D0->Kpi channel
 706
 707   //mass histograms
 708   Int_t pdgdaughtersD0[2]={211,321};//pi,K
 709   Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
 710   Double_t masses[2];
 711   masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
 712   masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
 713
 714   if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) fMassHist[index]->Fill(masses[0]);
 715   if(isSel>=2 && fPartOrAndAntiPart<=0) fMassHist[index]->Fill(masses[1]);
 716
 717
 718
 719   //MC histograms
 720   if(fReadMC){
 721
 722     Int_t matchtoMC=-1;
 723
 724     //D0
 725     Int_t pdgdaughters[2];
 726     pdgdaughters[0]=211;//pi
 727     pdgdaughters[1]=321;//K
 728
 729     matchtoMC = d->MatchToMC(fPDGmother,arrayMC,fNProngs,pdgdaughters);
 730
 731     Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
 732     if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0){ //D0
 733       if(matchtoMC>=0){
 734         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 735         Int_t pdgMC = dMC->GetPdgCode();
 736
 737         if(pdgMC==prongPdgPlus) fSigHist[index]->Fill(masses[0]);
 738         else fRflHist[index]->Fill(masses[0]);
 739
 740       } else fBkgHist[index]->Fill(masses[0]);
 741
 742     }
 743     if(isSel>=2 && fPartOrAndAntiPart<=0){ //D0bar
 744       if(matchtoMC>=0){
 745         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 746         Int_t pdgMC = dMC->GetPdgCode();
 747
 748         if(pdgMC==prongPdgMinus) fSigHist[index]->Fill(masses[1]);
 749         else fRflHist[index]->Fill(masses[1]);
 750       } else fBkgHist[index]->Fill(masses[1]);
 751     }
 752   }
 753 }
 754
 755 void AliAnalysisTaskSESignificance::FillDstar(AliAODRecoCascadeHF* dstarD0pi,TClonesArray *arrayMC,Int_t index,Int_t isSel){
 756     //D* channel
 757
 758     AliInfo("Dstar selected\n");
 759
 760     Double_t mass = dstarD0pi->DeltaInvMass();
 761
 762     if((isSel>0) && TMath::Abs(fPartOrAndAntiPart)>=0) fMassHist[index]->Fill(mass);
 763
 764     if(fReadMC) {
 765        Int_t matchtoMC = -1;
 766        matchtoMC = dstarD0pi->MatchToMC(413,421,fPDGDStarToD0pi, fPDGD0ToKpi, arrayMC);
 767
 768        Int_t prongPdgDStarPlus=413;//,prongPdgDStarMinus=(-1)*413;
 769
 770        if ((isSel>1) && TMath::Abs(fPartOrAndAntiPart)>=0) {
 771           //D*+
 772           if(matchtoMC>=0) {
 773              AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 774              Int_t pdgMC = dMC->GetPdgCode();
 775
 776              if (pdgMC==prongPdgDStarPlus) fSigHist[index]->Fill(mass);
 777              else {
 778                 dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
 779                 mass =  dstarD0pi->DeltaInvMass();
 780                 fRflHist[index]->Fill(mass);
 781                 dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
 782               }
 783             }
 784           else fBkgHist[index]->Fill(mass);
 785         }
 786     }
 787 }
 788
 789
 790 void AliAnalysisTaskSESignificance::FillDs(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel,Int_t optDecay){
 791
 792   // Fill Ds histos
 793
 794
 795   Int_t pdgDsKKpi[3]={321,321,211};//K,K,pi
 796   Int_t pdgDspiKK[3]={211,321,321};//pi,K,K
 797   Double_t masses[2];
 798   masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgDsKKpi); //Ds
 799   masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgDspiKK); //Dsbar
 800
 801   Int_t labDs=-1;
 802   if(fReadMC){
 803     labDs = d->MatchToMC(431,arrayMC,3,pdgDsKKpi);
 804   }
 805
 806   Int_t isKKpi=isSel&1;
 807   Int_t ispiKK=isSel&2;
 808   Int_t isPhiKKpi=isSel&4;
 809   Int_t isPhipiKK=isSel&8;
 810   Int_t isK0starKKpi=isSel&16;
 811   Int_t isK0starpiKK=isSel&32;
 812
 813
 814   if(fDsChannel==kPhi && (isPhiKKpi==0 && isPhipiKK==0)) return;
 815   if(fDsChannel==kK0star && (isK0starKKpi==0 && isK0starpiKK==0)) return;
 816
 817   if (optDecay==1){
 818     if(isKKpi  && fPartOrAndAntiPart*d->GetCharge()>=0) {
 819       if(fDsChannel==kPhi && isPhiKKpi==0) return;
 820       if(fDsChannel==kK0star && isK0starKKpi==0) return;
 821
 822       fMassHist[index]->Fill(masses[0]);
 823
 824       if(fReadMC){
 825         if(labDs>=0){
 826           Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 827           AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 828           Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 829           if(pdgCode0==321) {
 830             fSigHist[index]->Fill(masses[0]); //signal
 831           }else{
 832             fRflHist[index]->Fill(masses[0]); //Reflected signal
 833           }
 834         }else{
 835           fBkgHist[index]->Fill(masses[0]); // Background
 836         }
 837       }
 838     }
 839   }
 840
 841   if (optDecay==0){
 842     if(ispiKK && fPartOrAndAntiPart*d->GetCharge()>=0){
 843       if(fDsChannel==kPhi && isPhipiKK==0) return;
 844       if(fDsChannel==kK0star && isK0starpiKK==0) return;
 845
 846       fMassHist[index]->Fill(masses[1]);
 847
 848       if(fReadMC){
 849         if(labDs>=0){
 850           Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 851           AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 852           Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 853           if(pdgCode0==211) {
 854             fSigHist[index]->Fill(masses[1]);
 855           }else{
 856             fRflHist[index]->Fill(masses[1]);
 857           }
 858         }else{
 859           fBkgHist[index]->Fill(masses[1]);
 860         }
 861       }
 862     }
 863   }
 864 }
 865
 866 void AliAnalysisTaskSESignificance::FillD04p(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Int_t /*matchtoMC*/){
 867   //D0->Kpipipi channel
 868   AliInfo("D0 in 4 prongs channel not implemented\n");
 869
 870 }
 871
 872 void AliAnalysisTaskSESignificance::FillLambdac(AliAODRecoDecayHF* d,TClonesArray* arrayMC,Int_t index,Int_t isSel){
 873
 874   // Mass hypothesis
 875   Double_t masses[2];
 876   Int_t pdgdaughtersLc[3]={2212,321,211}; //p,K,pi
 877   masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc);
 878   Int_t pdgdaughtersLc2[3]={211,321,2212}; //pi,K,p
 879   masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc2);
 880
 881
 882   if(fPartOrAndAntiPart==0 || fPartOrAndAntiPart==d->GetCharge()) {
 883
 884     // isSel=1 : p K pi ; isSel=2 : pi K p ;
 885     if(isSel==1 || isSel==3) fMassHist[index]->Fill(masses[0]);
 886     if(isSel>=2) fMassHist[index]->Fill(masses[1]);
 887
 888     // Check the MC truth
 889     if(fReadMC){
 890
 891       Int_t ispKpi = 0;
 892       Int_t ispiKp = 0;
 893       ispKpi = isSel&1;
 894       ispiKp = isSel&2;
 895       Int_t matchtoMC = -1;
 896       //
 897       Int_t pPDG = 2212; // p
 898       Int_t kPDG = 321;  // K
 899       Int_t piPDG = 211; // pi
 900       Int_t absPdgMom = 4122;
 901       matchtoMC = d->MatchToMC(absPdgMom,arrayMC,fNProngs,pdgdaughtersLc);
 902
 903       if(matchtoMC>=0){
 904
 905         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 906         Int_t pdgMC = dMC->GetPdgCode();
 907         if (TMath::Abs(pdgMC)!=absPdgMom) AliInfo("What's up, isn't it a lambdac ?!");
 908         Int_t labDau0 = ((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 909         AliAODMCParticle* p0 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 910         Int_t pdgCode0 = TMath::Abs(p0->GetPdgCode());
 911         Int_t labDau1 = ((AliAODTrack*)d->GetDaughter(1))->GetLabel();
 912         AliAODMCParticle* p1 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau1);
 913         Int_t pdgCode1 = TMath::Abs(p1->GetPdgCode());
 914         Int_t labDau2 = ((AliAODTrack*)d->GetDaughter(2))->GetLabel();
 915         AliAODMCParticle* p2 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau2);
 916         Int_t pdgCode2 = TMath::Abs(p2->GetPdgCode());
 917
 918         // Fill in the histograms in case of p K pi decays
 919         if(ispKpi==1){
 920           if(pdgCode0==pPDG && pdgCode1==kPDG && pdgCode2==piPDG){
 921             fSigHist[index]->Fill(masses[0]);
 922           } else {
 923             fRflHist[index]->Fill(masses[0]);
 924           }
 925         }
 926         // Fill in the histograms in case of pi K p decays
 927         if(ispiKp==2){
 928           if(pdgCode0==piPDG && pdgCode1==kPDG && pdgCode2==pPDG){
 929             fSigHist[index]->Fill(masses[1]);
 930           } else {
 931             fRflHist[index]->Fill(masses[1]);
 932           }
 933         }
 934       } else {
 935         if(ispKpi==1) fBkgHist[index]->Fill(masses[0]);
 936         if(ispiKp==2) fBkgHist[index]->Fill(masses[1]);
 937       }
 938     }
 939   }
 940
 941
 942 }
 943
 944
 945 //________________________________________________________________________
 946 void AliAnalysisTaskSESignificance::Terminate(Option_t */*option*/)
 947 {
 948   // Terminate analysis
 949   //
 950   if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
 951
 952
 953   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 954   if (!fOutput) {
 955     printf("ERROR: fOutput not available\n");
 956     return;
 957   }
 958
 959   fCutList = dynamic_cast<TList*> (GetOutputData(2));
 960   if (!fCutList) {
 961     printf("ERROR: fCutList not available\n");
 962     return;
 963   }
 964
 965   AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
 966   if (!mdvtmp){
 967     cout<<"multidimvec not found in TList"<<endl;
 968     fCutList->ls();
 969     return;
 970   }
 971   Int_t nHist=mdvtmp->GetNTotCells();
 972   TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
 973   Bool_t drawn=kFALSE;
 974   for(Int_t i=0;i<nHist;i++){
 975
 976     TString hisname;
 977     TString signame;
 978     TString bkgname;
 979     TString rflname;
 980
 981     hisname.Form("hMass_%d",i);
 982     signame.Form("hSig_%d",i);
 983     bkgname.Form("hBkg_%d",i);
 984     rflname.Form("hRfl_%d",i);
 985
 986     fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 987
 988     if (!drawn && fMassHist[i]->GetEntries() > 0){
 989       c1->cd();
 990       fMassHist[i]->Draw();
 991       drawn=kTRUE;
 992     }
 993
 994     if(fReadMC){
 995       fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
 996       fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
 997       if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi) fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(rflname.Data()));
 998     }
 999
1000   }
1001
1002
1003
1004
1005   return;
1006 }
1007 //_________________________________________________________________________________________________
1008 Int_t AliAnalysisTaskSESignificance::CheckOrigin(const AliAODMCParticle* mcPart, const TClonesArray* mcArray)const{
1009
1010         //
1011         // checking whether the very mother of the D0 is a charm or a bottom quark
1012         //
1013
1014         Int_t pdgGranma = 0;
1015         Int_t mother = 0;
1016         mother = mcPart->GetMother();
1017         Int_t istep = 0;
1018         while (mother >0 ){
1019                 istep++;
1020                 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1021                 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1022                 if(!mcGranma) break;
1023                 pdgGranma = mcGranma->GetPdgCode();
1024                 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1025                 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1026                 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1027                         break;
1028                 }
1029                 mother = mcGranma->GetMother();
1030         }
1031         return pdgGranma;
1032 }
1033 //_________________________________________________________________________________________________