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