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   fHistNEvents->GetXaxis()->SetBinLabel(7,"MC Cand from c");
 348   fHistNEvents->GetXaxis()->SetBinLabel(8,"MC Cand from b");
 349   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 350   fOutput->Add(fHistNEvents);
 351
 352
 353   return;
 354 }
 355
 356 //________________________________________________________________________
 357 void AliAnalysisTaskSESignificance::UserExec(Option_t */*option*/)
 358 {
 359   // Execute analysis for current event:
 360   // heavy flavor candidates association to MC truth
 361
 362   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 363   if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
 364   // Post the data already here
 365   PostData(1,fOutput);
 366   TClonesArray *arrayProng =0;
 367
 368   if(!aod && AODEvent() && IsStandardAOD()) {
 369     // In case there is an AOD handler writing a standard AOD, use the AOD
 370     // event in memory rather than the input (ESD) event.
 371     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 372     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 373     // have to taken from the AOD event hold by the AliAODExtension
 374     AliAODHandler* aodHandler = (AliAODHandler*)
 375       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 376     if(aodHandler->GetExtensions()) {
 377
 378       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 379       AliAODEvent *aodFromExt = ext->GetAOD();
 380
 381
 382       switch(fDecChannel){
 383       case 0:
 384         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 385         break;
 386       case 1:
 387         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 388         break;
 389       case 2:
 390         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 391         break;
 392       case 3:
 393         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 394         break;
 395       case 4:
 396         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 397         break;
 398       case 5:
 399         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 400         break;
 401       }
 402     }
 403   } else {
 404     switch(fDecChannel){
 405     case 0:
 406       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 407       break;
 408     case 1:
 409       arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 410       break;
 411     case 2:
 412       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
 413       break;
 414     case 3:
 415       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 416       break;
 417     case 4:
 418       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm4Prong");
 419       break;
 420     case 5:
 421       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 422       break;
 423     }
 424   }
 425   if(!aod || !arrayProng) {
 426     AliError("AliAnalysisTaskSESignificance::UserExec:Branch not found!\n");
 427     return;
 428   }
 429
 430   // fix for temporary bug in ESDfilter
 431   // the AODs with null vertex pointer didn't pass the PhysSel
 432  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 433   TClonesArray *arrayMC=0;
 434   AliAODMCHeader *mcHeader=0;
 435
 436   // load MC particles
 437   if(fReadMC){
 438
 439     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 440     if(!arrayMC) {
 441       AliWarning("AliAnalysisTaskSESignificance::UserExec:MC particles branch not found!\n");
 442       //    return;
 443     }
 444
 445     // load MC header
 446     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 447     if(!mcHeader) {
 448       AliError("AliAnalysisTaskSESignificance::UserExec:MC header branch not found!\n");
 449       return;
 450     }
 451   }
 452
 453   fHistNEvents->Fill(0); // count event
 454
 455   AliAODRecoDecayHF *d=0;
 456   Int_t nprongs=1;
 457   Int_t *pdgdaughters=0x0;
 458   Int_t absPdgMom=411;
 459
 460
 461   switch(fDecChannel){
 462   case 0:
 463     //D+
 464     pdgdaughters =new Int_t[3];
 465     pdgdaughters[0]=211;//pi
 466     pdgdaughters[1]=321;//K
 467     pdgdaughters[2]=211;//pi
 468     nprongs=3;
 469     absPdgMom=411;
 470     break;
 471   case 1:
 472     //D0
 473     pdgdaughters =new Int_t[2];
 474     pdgdaughters[0]=211;//pi
 475     pdgdaughters[1]=321;//K
 476     nprongs=2;
 477     absPdgMom=421;
 478     break;
 479   case 2:
 480     //D*
 481     pdgdaughters =new Int_t[3];
 482     pdgdaughters[1]=211;//pi
 483     pdgdaughters[0]=321;//K
 484     pdgdaughters[2]=211;//pi (soft?)
 485     nprongs=3;
 486     absPdgMom=413;
 487     break;
 488   case 3:
 489     //Ds
 490     pdgdaughters =new Int_t[3];
 491     pdgdaughters[0]=321;//K
 492     pdgdaughters[1]=321;//K
 493     pdgdaughters[2]=211;//pi
 494     nprongs=3;
 495     absPdgMom=431;
 496     break;
 497   case 4:
 498     //D0 in 4 prongs
 499     pdgdaughters =new Int_t[4];
 500     pdgdaughters[0]=321;
 501     pdgdaughters[1]=211;
 502     pdgdaughters[2]=211;
 503     pdgdaughters[3]=211;
 504     nprongs=4;
 505     absPdgMom=421;
 506     break;
 507   case 5:
 508     //Lambda_c
 509     pdgdaughters =new Int_t[3];
 510     pdgdaughters[0]=2212;//p
 511     pdgdaughters[1]=321;//K
 512     pdgdaughters[2]=211;//pi
 513     nprongs=3;
 514     absPdgMom=4122;
 515     break;
 516   }
 517
 518   Int_t nHistpermv=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
 519   Int_t nProng = arrayProng->GetEntriesFast();
 520   if(fDebug>1) printf("Number of D2H: %d\n",nProng);
 521
 522   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
 523   TString trigclass=aod->GetFiredTriggerClasses();
 524   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(5);
 525
 526   if(fRDCuts->IsEventSelected(aod)) fHistNEvents->Fill(1);
 527   else{
 528     if(fRDCuts->GetWhyRejection()==1) // rejected for pileup
 529       fHistNEvents->Fill(4);
 530     return;
 531   }
 532
 533   for (Int_t iProng = 0; iProng < nProng; iProng++) {
 534
 535     d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
 536
 537     Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(absPdgMom));
 538     Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel,aod);
 539
 540     if(fReadMC && fBFeedDown!=kBoth && isSelected){
 541       Int_t labD = d->MatchToMC(absPdgMom,arrayMC,nprongs,pdgdaughters);
 542       if(labD>=0){
 543         AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
 544         Int_t label=partD->GetMother();
 545         AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(label);
 546         while(label>=0){//get first mother
 547           mot = (AliAODMCParticle*)arrayMC->At(label);
 548           label=mot->GetMother();
 549         }
 550         Int_t pdgMotCode = mot->GetPdgCode();
 551
 552         if(TMath::Abs(pdgMotCode)<=4){
 553           fHistNEvents->Fill(6);
 554           if(fBFeedDown==kBeautyOnly)isSelected=kFALSE; //from primary charm
 555         }else{
 556           fHistNEvents->Fill(7);
 557           if(fBFeedDown==kCharmOnly) isSelected=kFALSE; //from beauty
 558         }
 559
 560         /*
 561         if(TMath::Abs(pdgMotCode)==4 && fBFeedDown==kBeautyOnly) isSelected=kFALSE; //from primary charm
 562         if(TMath::Abs(pdgMotCode)==5 && fBFeedDown==kCharmOnly) isSelected=kFALSE; //from beauty
 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     }// end if selected
 621
 622   }
 623
 624   if(pdgdaughters) {delete [] pdgdaughters; pdgdaughters=NULL;}
 625
 626   PostData(1,fOutput);
 627   return;
 628 }
 629
 630 //***************************************************************************
 631
 632 // Methods used in the UserExec
 633
 634 void AliAnalysisTaskSESignificance::CalculateInvMasses(AliAODRecoDecayHF* d,Double_t*& masses,Int_t& nmasses){
 635   //Calculates all the possible invariant masses for each candidate
 636   //NB: the implementation for each candidate is responsibility of the corresponding developer
 637
 638   switch(fDecChannel){
 639   case 0:
 640     //D+ -- Giacomo, Renu ?
 641     {
 642       nmasses=1;
 643       masses=new Double_t[nmasses];
 644       Int_t pdgdaughters[3] = {211,321,211};
 645       masses[0]=d->InvMass(3,(UInt_t*)pdgdaughters);
 646     }
 647     break;
 648   case 1:
 649     //D0 (Kpi)  -- Chiara
 650     {
 651       const Int_t ndght=2;
 652       nmasses=2;
 653       masses=new Double_t[nmasses];
 654       Int_t pdgdaughtersD0[ndght]={211,321};//pi,K
 655       masses[0]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0); //D0
 656       Int_t pdgdaughtersD0bar[ndght]={321,211};//K,pi
 657       masses[1]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0bar); //D0bar
 658     }
 659     break;
 660   case 2:
 661     //D* -- ? Yifei, Alessandro
 662     {
 663       //.....
 664     }
 665     break;
 666   case 3:
 667     //Ds  -- Sergey, Sahdana
 668     {
 669       const Int_t ndght=3;
 670       nmasses=2;
 671       masses=new Double_t[nmasses];
 672       Int_t pdgdaughtersDsKKpi[ndght]={321,321,211};//K,K,pi
 673       masses[0]=d->InvMass(ndght,(UInt_t*)pdgdaughtersDsKKpi); //Ds
 674       Int_t pdgdaughtersDspiKK[ndght]={211,321,321};//pi,K,K
 675       masses[1]=d->InvMass(ndght,(UInt_t*)pdgdaughtersDspiKK); //Dsbar
 676
 677       //.....
 678     }
 679     break;
 680   case 4:
 681     //D0 (Kpipipi) -- ? Rossella, Fabio ???
 682     {
 683       //.....
 684     }
 685     break;
 686   case 5:
 687     //Lambda_c -- Rossella
 688     {
 689       //.....
 690     }
 691     break;
 692   default:
 693     break;
 694   }
 695 }
 696
 697 //********************************************************************************************
 698
 699 //Methods to fill the istograms with MC information, one for each candidate
 700 //NB: the implementation for each candidate is responsibility of the corresponding developer
 701
 702 void AliAnalysisTaskSESignificance::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Double_t* masses,Int_t isSel){
 703     //D+ channel
 704   if(!isSel){
 705     AliError("Candidate not selected\n");
 706     return;
 707   }
 708
 709   if(fPartOrAndAntiPart*d->GetCharge()<0)return;
 710
 711   fMassHist[index]->Fill(masses[0]);
 712
 713   Int_t pdgdaughters[3] = {211,321,211};
 714
 715   if(fReadMC){
 716     Int_t lab=-1;
 717     lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
 718     if(lab>=0){ //signal
 719       fSigHist[index]->Fill(masses[0]);
 720     } else{ //background
 721       fBkgHist[index]->Fill(masses[0]);
 722     }
 723   }
 724 }
 725
 726
 727 void AliAnalysisTaskSESignificance::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Double_t* masses,Int_t isSel){
 728
 729   //D0->Kpi channel
 730
 731   //mass histograms
 732   if(!masses){
 733     AliError("Masses not calculated\n");
 734     return;
 735   }
 736   if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) fMassHist[index]->Fill(masses[0]);
 737   if(isSel>=2 && fPartOrAndAntiPart<=0) fMassHist[index]->Fill(masses[1]);
 738
 739
 740
 741   //MC histograms
 742   if(fReadMC){
 743
 744     Int_t matchtoMC=-1;
 745
 746     //D0
 747     Int_t pdgdaughters[2];
 748     pdgdaughters[0]=211;//pi
 749     pdgdaughters[1]=321;//K
 750     Int_t nprongs=2;
 751     Int_t absPdgMom=421;
 752
 753     matchtoMC = d->MatchToMC(absPdgMom,arrayMC,nprongs,pdgdaughters);
 754
 755     Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
 756     if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0){ //D0
 757       if(matchtoMC>=0){
 758         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 759         Int_t pdgMC = dMC->GetPdgCode();
 760
 761         if(pdgMC==prongPdgPlus) fSigHist[index]->Fill(masses[0]);
 762         else fRflHist[index]->Fill(masses[0]);
 763
 764       } else fBkgHist[index]->Fill(masses[0]);
 765
 766     }
 767     if(isSel>=2 && fPartOrAndAntiPart<=0){ //D0bar
 768       if(matchtoMC>=0){
 769         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 770         Int_t pdgMC = dMC->GetPdgCode();
 771
 772         if(pdgMC==prongPdgMinus) fSigHist[index]->Fill(masses[1]);
 773         else fRflHist[index]->Fill(masses[1]);
 774       } else fBkgHist[index]->Fill(masses[1]);
 775     }
 776   }
 777 }
 778
 779 void AliAnalysisTaskSESignificance::FillDstar(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Double_t* /*masses*/,Int_t /*matchtoMC*/){
 780     //D* channel
 781
 782   AliInfo("Dstar channel not implemented\n");
 783
 784 }
 785
 786
 787 void AliAnalysisTaskSESignificance::FillDs(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Double_t* masses,Int_t isSel){
 788
 789   //AliInfo("Ds channel not implemented\n");
 790
 791
 792   if(!masses){
 793     AliError("Masses not calculated\n");
 794     return;
 795   }
 796
 797   Int_t pdgDstoKKpi[3]={321,321,211};
 798   Int_t labDs=-1;
 799   if(fReadMC){
 800     labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
 801   }
 802
 803
 804   if(isSel&1  && fPartOrAndAntiPart*d->GetCharge()>=0) {
 805
 806     fMassHist[index]->Fill(masses[0]);
 807
 808     if(fReadMC){
 809
 810       if(labDs>=0){
 811         Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 812         AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 813         Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 814
 815         if(pdgCode0==321) {
 816
 817           fSigHist[index]->Fill(masses[0]); //signal
 818         }else{
 819           fRflHist[index]->Fill(masses[0]); //Reflected signal
 820         }
 821       }else{
 822         fBkgHist[index]->Fill(masses[0]); // Background
 823       }
 824     }
 825   }
 826
 827   if(isSel&2 && fPartOrAndAntiPart*d->GetCharge()>=0){
 828     fMassHist[index]->Fill(masses[1]);
 829     if(fReadMC){
 830       if(labDs>=0){
 831         Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 832         AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 833         Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 834
 835         if(pdgCode0==211) {
 836
 837           fSigHist[index]->Fill(masses[1]);
 838         }else{
 839           fRflHist[index]->Fill(masses[1]);
 840         }
 841       }else{
 842         fBkgHist[index]->Fill(masses[1]);
 843       }
 844     }
 845   }
 846
 847
 848   //MC histograms
 849
 850   //Ds channel
 851         // Int_t isKKpi=0;
 852         // Int_t ispiKK=0;
 853         //   isKKpi=isSelected&1;
 854         //   ispiKK=isSelected&2;
 855
 856         // Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 857         // AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 858         // Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 859
 860         //       if(isKKpi){
 861         //        if(pdgCode0==321){
 862         //         fSigHist[index]->Fill(masses[0]);
 863         //        }else{
 864         //         fRflHist[index]->Fill(masses[0]);
 865         //        }
 866         //       }
 867         //       if(ispiKK){
 868         //        if(pdgCode0==211){
 869         //         fSigHist[index]->Fill(masses[1]);
 870         //        }else{
 871         //         fRflHist[index]->Fill(masses[1]);
 872         //        }
 873         //       }
 874
 875 }
 876
 877 void AliAnalysisTaskSESignificance::FillD04p(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Double_t* /*masses*/,Int_t /*matchtoMC*/){
 878   //D0->Kpipipi channel
 879   AliInfo("D0 in 4 prongs channel not implemented\n");
 880
 881 }
 882
 883 void AliAnalysisTaskSESignificance::FillLambdac(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Double_t* /*masses*/,Int_t /*matchtoMC*/){
 884   AliInfo("Lambdac channel not implemented\n");
 885
 886   //Lambdac channel
 887   // Int_t ispKpi=0;
 888   // Int_t ispiKp=0;
 889
 890   // ispKpi=isSelected&1;
 891   // ispiKp=isSelected&2;
 892   // if(matchtoMC>=0){
 893   // Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 894   // AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 895   // Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 896   // if(ispKpi){
 897   //   if(pdgCode0==2212){
 898   //     fSigHist[index]->Fill(invMass[0]);
 899   //   }else{
 900   //     fRflHist[index]->Fill(invMass[0]);
 901   //   }
 902   // }
 903   // if(ispiKp){
 904   //   if(pdgCode0==211){
 905   //     fSigHist[index]->Fill(invMass[1]);
 906   //   }else{
 907   //     fRflHist[index]->Fill(invMass[1]);
 908   //   }
 909   // }
 910   // }else{
 911   //   fBkgHist[index]
 912   // }
 913
 914
 915 }
 916
 917
 918 //________________________________________________________________________
 919 void AliAnalysisTaskSESignificance::Terminate(Option_t */*option*/)
 920 {
 921   // Terminate analysis
 922   //
 923   if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
 924
 925
 926   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 927   if (!fOutput) {
 928     printf("ERROR: fOutput not available\n");
 929     return;
 930   }
 931
 932   fCutList = dynamic_cast<TList*> (GetOutputData(2));
 933   if (!fCutList) {
 934     printf("ERROR: fCutList not available\n");
 935     return;
 936   }
 937
 938   AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
 939   if (!mdvtmp){
 940     cout<<"multidimvec not found in TList"<<endl;
 941     fCutList->ls();
 942     return;
 943   }
 944   Int_t nHist=mdvtmp->GetNTotCells();
 945   TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
 946   Bool_t drawn=kFALSE;
 947   for(Int_t i=0;i<nHist;i++){
 948
 949     TString hisname;
 950     TString signame;
 951     TString bkgname;
 952     TString rflname;
 953
 954     hisname.Form("hMass_%d",i);
 955     signame.Form("hSig_%d",i);
 956     bkgname.Form("hBkg_%d",i);
 957     rflname.Form("hRfl_%d",i);
 958
 959     fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 960
 961     if (!drawn && fMassHist[i]->GetEntries() > 0){
 962       c1->cd();
 963       fMassHist[i]->Draw();
 964       drawn=kTRUE;
 965     }
 966
 967     if(fReadMC){
 968       fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
 969       fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
 970       if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi) fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(rflname.Data()));
 971     }
 972
 973   }
 974
 975
 976
 977
 978   return;
 979 }
 980 //-------------------------------------------
 981