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