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==1 || isSel==3) && fPartOrAndAntiPart>=0) fMassHist[index]->Fill(mass);
 763     if(isSel>=2 && fPartOrAndAntiPart<=0) fMassHist[index]->Fill(mass);
 764
 765     if(fReadMC) {
 766        Int_t matchtoMC = -1;
 767        matchtoMC = dstarD0pi->MatchToMC(413,421,fPDGDStarToD0pi, fPDGD0ToKpi, arrayMC);
 768
 769        Int_t prongPdgDStarPlus=413,prongPdgDStarMinus=(-1)*413;
 770
 771        if ((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) {
 772           //D*+
 773           if(matchtoMC>=0) {
 774              AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 775              Int_t pdgMC = dMC->GetPdgCode();
 776
 777              if (pdgMC==prongPdgDStarPlus) fSigHist[index]->Fill(mass);
 778              else {
 779                 dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
 780                 mass =  dstarD0pi->DeltaInvMass();
 781                 fRflHist[index]->Fill(mass);
 782                 dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
 783               }
 784             }
 785           else fBkgHist[index]->Fill(mass);
 786         }
 787
 788         if (isSel>=2 && fPartOrAndAntiPart<=0) {
 789            //D*-
 790            if (matchtoMC>=0) {
 791               AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 792               Int_t pdgMC = dMC->GetPdgCode();
 793
 794               if (pdgMC==prongPdgDStarMinus) fSigHist[index]->Fill(mass);
 795               else {
 796                  dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
 797                  mass = dstarD0pi->DeltaInvMass();
 798                  fRflHist[index]->Fill(mass);
 799                  dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
 800               }
 801             }
 802             else fBkgHist[index]->Fill(mass);
 803           }
 804        }
 805
 806 }
 807
 808
 809 void AliAnalysisTaskSESignificance::FillDs(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel,Int_t optDecay){
 810
 811   // Fill Ds histos
 812
 813
 814   Int_t pdgDsKKpi[3]={321,321,211};//K,K,pi
 815   Int_t pdgDspiKK[3]={211,321,321};//pi,K,K
 816   Double_t masses[2];
 817   masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgDsKKpi); //Ds
 818   masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgDspiKK); //Dsbar
 819
 820   Int_t labDs=-1;
 821   if(fReadMC){
 822     labDs = d->MatchToMC(431,arrayMC,3,pdgDsKKpi);
 823   }
 824
 825   Int_t isKKpi=isSel&1;
 826   Int_t ispiKK=isSel&2;
 827   Int_t isPhiKKpi=isSel&4;
 828   Int_t isPhipiKK=isSel&8;
 829   Int_t isK0starKKpi=isSel&16;
 830   Int_t isK0starpiKK=isSel&32;
 831
 832
 833   if(fDsChannel==kPhi && (isPhiKKpi==0 && isPhipiKK==0)) return;
 834   if(fDsChannel==kK0star && (isK0starKKpi==0 && isK0starpiKK==0)) return;
 835
 836   if (optDecay==1){
 837     if(isKKpi  && fPartOrAndAntiPart*d->GetCharge()>=0) {
 838       if(fDsChannel==kPhi && isPhiKKpi==0) return;
 839       if(fDsChannel==kK0star && isK0starKKpi==0) return;
 840
 841       fMassHist[index]->Fill(masses[0]);
 842
 843       if(fReadMC){
 844         if(labDs>=0){
 845           Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 846           AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 847           Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 848           if(pdgCode0==321) {
 849             fSigHist[index]->Fill(masses[0]); //signal
 850           }else{
 851             fRflHist[index]->Fill(masses[0]); //Reflected signal
 852           }
 853         }else{
 854           fBkgHist[index]->Fill(masses[0]); // Background
 855         }
 856       }
 857     }
 858   }
 859
 860   if (optDecay==0){
 861     if(ispiKK && fPartOrAndAntiPart*d->GetCharge()>=0){
 862       if(fDsChannel==kPhi && isPhipiKK==0) return;
 863       if(fDsChannel==kK0star && isK0starpiKK==0) return;
 864
 865       fMassHist[index]->Fill(masses[1]);
 866
 867       if(fReadMC){
 868         if(labDs>=0){
 869           Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 870           AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 871           Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 872           if(pdgCode0==211) {
 873             fSigHist[index]->Fill(masses[1]);
 874           }else{
 875             fRflHist[index]->Fill(masses[1]);
 876           }
 877         }else{
 878           fBkgHist[index]->Fill(masses[1]);
 879         }
 880       }
 881     }
 882   }
 883 }
 884
 885 void AliAnalysisTaskSESignificance::FillD04p(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Int_t /*matchtoMC*/){
 886   //D0->Kpipipi channel
 887   AliInfo("D0 in 4 prongs channel not implemented\n");
 888
 889 }
 890
 891 void AliAnalysisTaskSESignificance::FillLambdac(AliAODRecoDecayHF* d,TClonesArray* arrayMC,Int_t index,Int_t isSel){
 892
 893   // Mass hypothesis
 894   Double_t masses[2];
 895   Int_t pdgdaughtersLc[3]={2212,321,211}; //p,K,pi
 896   masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc);
 897   Int_t pdgdaughtersLc2[3]={211,321,2212}; //pi,K,p
 898   masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc2);
 899
 900
 901   if(fPartOrAndAntiPart==0 || fPartOrAndAntiPart==d->GetCharge()) {
 902
 903     // isSel=1 : p K pi ; isSel=2 : pi K p ;
 904     if(isSel==1 || isSel==3) fMassHist[index]->Fill(masses[0]);
 905     if(isSel>=2) fMassHist[index]->Fill(masses[1]);
 906
 907     // Check the MC truth
 908     if(fReadMC){
 909
 910       Int_t ispKpi = 0;
 911       Int_t ispiKp = 0;
 912       ispKpi = isSel&1;
 913       ispiKp = isSel&2;
 914       Int_t matchtoMC = -1;
 915       //
 916       Int_t pPDG = 2212; // p
 917       Int_t kPDG = 321;  // K
 918       Int_t piPDG = 211; // pi
 919       Int_t absPdgMom = 4122;
 920       matchtoMC = d->MatchToMC(absPdgMom,arrayMC,fNProngs,pdgdaughtersLc);
 921
 922       if(matchtoMC>=0){
 923
 924         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 925         Int_t pdgMC = dMC->GetPdgCode();
 926         if (TMath::Abs(pdgMC)!=absPdgMom) AliInfo("What's up, isn't it a lambdac ?!");
 927         Int_t labDau0 = ((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 928         AliAODMCParticle* p0 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 929         Int_t pdgCode0 = TMath::Abs(p0->GetPdgCode());
 930         Int_t labDau1 = ((AliAODTrack*)d->GetDaughter(1))->GetLabel();
 931         AliAODMCParticle* p1 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau1);
 932         Int_t pdgCode1 = TMath::Abs(p1->GetPdgCode());
 933         Int_t labDau2 = ((AliAODTrack*)d->GetDaughter(2))->GetLabel();
 934         AliAODMCParticle* p2 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau2);
 935         Int_t pdgCode2 = TMath::Abs(p2->GetPdgCode());
 936
 937         // Fill in the histograms in case of p K pi decays
 938         if(ispKpi==1){
 939           if(pdgCode0==pPDG && pdgCode1==kPDG && pdgCode2==piPDG){
 940             fSigHist[index]->Fill(masses[0]);
 941           } else {
 942             fRflHist[index]->Fill(masses[0]);
 943           }
 944         }
 945         // Fill in the histograms in case of pi K p decays
 946         if(ispiKp==2){
 947           if(pdgCode0==piPDG && pdgCode1==kPDG && pdgCode2==pPDG){
 948             fSigHist[index]->Fill(masses[1]);
 949           } else {
 950             fRflHist[index]->Fill(masses[1]);
 951           }
 952         }
 953       } else {
 954         if(ispKpi==1) fBkgHist[index]->Fill(masses[0]);
 955         if(ispiKp==2) fBkgHist[index]->Fill(masses[1]);
 956       }
 957     }
 958   }
 959
 960
 961 }
 962
 963
 964 //________________________________________________________________________
 965 void AliAnalysisTaskSESignificance::Terminate(Option_t */*option*/)
 966 {
 967   // Terminate analysis
 968   //
 969   if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
 970
 971
 972   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 973   if (!fOutput) {
 974     printf("ERROR: fOutput not available\n");
 975     return;
 976   }
 977
 978   fCutList = dynamic_cast<TList*> (GetOutputData(2));
 979   if (!fCutList) {
 980     printf("ERROR: fCutList not available\n");
 981     return;
 982   }
 983
 984   AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
 985   if (!mdvtmp){
 986     cout<<"multidimvec not found in TList"<<endl;
 987     fCutList->ls();
 988     return;
 989   }
 990   Int_t nHist=mdvtmp->GetNTotCells();
 991   TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
 992   Bool_t drawn=kFALSE;
 993   for(Int_t i=0;i<nHist;i++){
 994
 995     TString hisname;
 996     TString signame;
 997     TString bkgname;
 998     TString rflname;
 999
1000     hisname.Form("hMass_%d",i);
1001     signame.Form("hSig_%d",i);
1002     bkgname.Form("hBkg_%d",i);
1003     rflname.Form("hRfl_%d",i);
1004
1005     fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1006
1007     if (!drawn && fMassHist[i]->GetEntries() > 0){
1008       c1->cd();
1009       fMassHist[i]->Draw();
1010       drawn=kTRUE;
1011     }
1012
1013     if(fReadMC){
1014       fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
1015       fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
1016       if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi) fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(rflname.Data()));
1017     }
1018
1019   }
1020
1021
1022
1023
1024   return;
1025 }
1026 //_________________________________________________________________________________________________
1027 Int_t AliAnalysisTaskSESignificance::CheckOrigin(const AliAODMCParticle* mcPart, const TClonesArray* mcArray)const{
1028
1029         //
1030         // checking whether the very mother of the D0 is a charm or a bottom quark
1031         //
1032
1033         Int_t pdgGranma = 0;
1034         Int_t mother = 0;
1035         mother = mcPart->GetMother();
1036         Int_t istep = 0;
1037         while (mother >0 ){
1038                 istep++;
1039                 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1040                 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1041                 if(!mcGranma) break;
1042                 pdgGranma = mcGranma->GetPdgCode();
1043                 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1044                 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1045                 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1046                         break;
1047                 }
1048                 mother = mcGranma->GetMother();
1049         }
1050         return pdgGranma;
1051 }
1052 //_________________________________________________________________________________________________