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