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