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