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