PWG3/vertexingHF/AliAnalysisTaskSESignificance.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /////////////////////////////////////////////////////////////
  17 //
  18 // AliAnalysisTaskSESignificane to calculate effects on
  19 // significance of D mesons  cut
  20 // Authors: G. Ortona, ortona@to.infn.it
  21 // F. Prino, prino@to.infn.it
  22 // Renu Bala, bala@to.infn.it
  23 // Chiara Bianchin, cbianchi@pd.infn.it
  24 /////////////////////////////////////////////////////////////
  25
  26 #include <Riostream.h>
  27 #include <TClonesArray.h>
  28 #include <TCanvas.h>
  29 #include <TList.h>
  30 #include <TFile.h>
  31 #include <TString.h>
  32 #include <TH1F.h>
  33 #include <TDatabasePDG.h>
  34
  35 #include <AliLog.h>
  36 #include "AliAnalysisManager.h"
  37 #include "AliAODHandler.h"
  38 #include "AliAODEvent.h"
  39 #include "AliAODVertex.h"
  40 #include "AliAODTrack.h"
  41 #include "AliAODMCHeader.h"
  42 #include "AliAODMCParticle.h"
  43 #include "AliAODRecoDecayHF3Prong.h"
  44 #include "AliAODRecoDecayHF.h"
  45 #include "AliAODRecoDecayHF2Prong.h"
  46 #include "AliAODRecoDecayHF4Prong.h"
  47 #include "AliAODRecoCascadeHF.h"
  48
  49 #include "AliAnalysisTaskSE.h"
  50 #include "AliRDHFCutsDplustoKpipi.h"
  51 #include "AliRDHFCutsD0toKpipipi.h"
  52 #include "AliRDHFCutsDstoKKpi.h"
  53 #include "AliRDHFCutsDStartoKpipi.h"
  54 #include "AliRDHFCutsD0toKpi.h"
  55 #include "AliRDHFCutsLctopKpi.h"
  56 #include "AliMultiDimVector.h"
  57
  58 #include "AliAnalysisTaskSESignificance.h"
  59
  60 ClassImp(AliAnalysisTaskSESignificance)
  61
  62
  63 //________________________________________________________________________
  64 AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance():
  65   AliAnalysisTaskSE(),
  66   fOutput(0),
  67   fCutList(0),
  68   fHistNEvents(0),
  69   fUpmasslimit(1.965),
  70   fLowmasslimit(1.765),
  71   fRDCuts(0),
  72   fNPtBins(0),
  73   fReadMC(kFALSE),
  74   fDecChannel(0),
  75   fSelectionlevel(0)
  76 {
  77   // Default constructor
  78 }
  79
  80 //________________________________________________________________________
  81 AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance(const char *name, TList* listMDV,AliRDHFCuts *rdCuts,Int_t decaychannel,Int_t selectionlevel):
  82   AliAnalysisTaskSE(name),
  83   fOutput(0),
  84   fCutList(listMDV),
  85   fHistNEvents(0),
  86   fUpmasslimit(0),
  87   fLowmasslimit(0),
  88   fRDCuts(rdCuts),
  89   fNPtBins(0),
  90   fReadMC(kFALSE),
  91   fDecChannel(decaychannel),
  92   fSelectionlevel(selectionlevel)
  93 {
  94
  95   Int_t pdg=421;
  96   switch(fDecChannel){
  97   case 0:
  98     pdg=411;
  99     break;
 100   case 1:
 101     pdg=421;
 102     break;
 103   case 2:
 104     pdg=413;
 105     break;
 106   case 3:
 107     pdg=431;
 108     break;
 109   case 4:
 110     pdg=421;
 111     break;
 112   case 5:
 113     pdg=4122;
 114     break;
 115   }
 116
 117   SetMassLimits(0.1,pdg); //check range
 118   fNPtBins=fRDCuts->GetNPtBins();
 119
 120   if(fDebug>1)fRDCuts->PrintAll();
 121    // Output slot #1 writes into a TList container
 122   DefineOutput(1,TList::Class());  //My private output
 123   DefineOutput(2,TList::Class());
 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 /*
 145   if(fRDCuts) {
 146     delete fRDCuts;
 147     fRDCuts= 0;
 148   }
 149 */
 150
 151 }
 152 //_________________________________________________________________
 153 Bool_t AliAnalysisTaskSESignificance::CheckConsistency(){
 154
 155   Bool_t result = kTRUE;
 156
 157   const Int_t nvars=fRDCuts->GetNVars();//ForOpt();
 158   //Float_t *vars = new Float_t[nvars];
 159   Bool_t *varsforopt = fRDCuts->GetVarsForOpt();
 160   Bool_t *uppervars = fRDCuts->GetIsUpperCut();
 161   TString *names = fRDCuts->GetVarNames();
 162
 163   for(Int_t i=0;i<fNPtBins;i++){
 164     TString mdvname=Form("multiDimVectorPtBin%d",i);
 165     Int_t ic=0;
 166
 167     for(Int_t ivar=0;ivar<nvars;ivar++){
 168       if(varsforopt[ivar]){
 169         Float_t min = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMinLimit(ic);
 170         Float_t max = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMaxLimit(ic);
 171         if(min==max){
 172           printf("AliAnalysisTaskSESignificance::CheckConsistency: ERROR! \n tight and loose cut for optimization variable number %d are the same in ptbin %d\n",ic,i);
 173           result = kFALSE;
 174         }
 175         Bool_t lowermdv = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGreaterThan(ic);
 176         if(uppervars[ivar]&&lowermdv){
 177           printf("AliAnalysisTaskSESignificance::CheckConsistency: WARNING! %s is declared as uppercut, but as been given tighter cut larger then loose cut in ptbin %d \n ---Task will use swapped Tight/Loose cuts \n ",names[ivar].Data(),i);
 178           ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->SwapLimits(ic);
 179           result = kTRUE;
 180         }
 181         if(!uppervars[ivar]&&!lowermdv){
 182           printf("AliAnalysisTaskSESignificance::CheckConsistency: WARNING! %s is declared as lower cut, but as been given tighter cut smaller then loose cut in ptbin %d \n ---Task will use swapped Tight/Loose cuts \n",names[ivar].Data(),i);
 183           ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->SwapLimits(ic);
 184           result = kTRUE;
 185         }
 186         ic++;
 187       }
 188     }
 189   }
 190   return result;
 191 }
 192 //_________________________________________________________________
 193 void  AliAnalysisTaskSESignificance::SetMassLimits(Float_t range, Int_t pdg){
 194   Float_t mass=0;
 195   Int_t abspdg=TMath::Abs(pdg);
 196   mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
 197   fUpmasslimit = mass+range;
 198   fLowmasslimit = mass-range;
 199 }
 200 //_________________________________________________________________
 201 void  AliAnalysisTaskSESignificance::SetMassLimits(Float_t lowlimit, Float_t uplimit){
 202   if(uplimit>lowlimit)
 203     {
 204       fUpmasslimit = lowlimit;
 205       fLowmasslimit = uplimit;
 206     }
 207 }
 208
 209
 210
 211 //________________________________________________________________________
 212 void AliAnalysisTaskSESignificance::LocalInit()
 213 {
 214   // Initialization
 215
 216   if(fDebug > 1) printf("AnalysisTaskSESignificance::Init() \n");
 217
 218   TList *mdvList =  new TList();
 219   mdvList->SetOwner();
 220   mdvList = fCutList;
 221
 222   PostData(2,mdvList);
 223
 224   return;
 225 }
 226 //________________________________________________________________________
 227 void AliAnalysisTaskSESignificance::UserCreateOutputObjects()
 228 {
 229   // Create the output container
 230
 231   if(fDebug > 1) printf("AnalysisTaskSESignificance::UserCreateOutputObjects() \n");
 232
 233   // Several histograms are more conveniently managed in a TList
 234   fOutput = new TList();
 235   fOutput->SetOwner();
 236   fOutput->SetName("OutputHistos");
 237
 238   //same number of steps in each multiDimVectorPtBin%d !
 239   Int_t nHist=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
 240   cout<<"ncells = "<<nHist<<" n ptbins = "<<fNPtBins<<endl;
 241   nHist=nHist*fNPtBins;
 242   cout<<"Total = "<<nHist<<endl;
 243   for(Int_t i=0;i<nHist;i++){
 244
 245     TString hisname;
 246     TString signame;
 247     TString bkgname;
 248     TString rflname;
 249     TString title;
 250
 251     hisname.Form("hMass_%d",i);
 252     signame.Form("hSig_%d",i);
 253     bkgname.Form("hBkg_%d",i);
 254     rflname.Form("hRfl_%d",i);
 255
 256     title.Form("Invariant mass;M[GeV/c^{2}];Entries");
 257
 258     fMassHist[i]=new TH1F(hisname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
 259     fMassHist[i]->Sumw2();
 260     fOutput->Add(fMassHist[i]);
 261
 262     if(fReadMC){
 263       fSigHist[i]=new TH1F(signame.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
 264       fSigHist[i]->Sumw2();
 265       fOutput->Add(fSigHist[i]);
 266
 267       fBkgHist[i]=new TH1F(bkgname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
 268       fBkgHist[i]->Sumw2();
 269       fOutput->Add(fBkgHist[i]);
 270
 271       fRflHist[i]=new TH1F(rflname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
 272       fRflHist[i]->Sumw2();
 273       fOutput->Add(fRflHist[i]);
 274
 275     }
 276   }
 277
 278   fHistNEvents=new TH1F("fHistNEvents","Number of AODs scanned",4,0,4.);
 279   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 280   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvNotSelected");
 281   fHistNEvents->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
 282   fHistNEvents->GetXaxis()->SetBinLabel(4,"nTotEntries Mass hists");
 283   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 284   fOutput->Add(fHistNEvents);
 285
 286
 287   return;
 288 }
 289
 290 //________________________________________________________________________
 291 void AliAnalysisTaskSESignificance::UserExec(Option_t */*option*/)
 292 {
 293   // Execute analysis for current event:
 294   // heavy flavor candidates association to MC truth
 295
 296   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 297   if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
 298   // Post the data already here
 299   PostData(1,fOutput);
 300   TClonesArray *arrayProng =0;
 301
 302   if(!aod && AODEvent() && IsStandardAOD()) {
 303     // In case there is an AOD handler writing a standard AOD, use the AOD
 304     // event in memory rather than the input (ESD) event.
 305     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 306     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 307     // have to taken from the AOD event hold by the AliAODExtension
 308     AliAODHandler* aodHandler = (AliAODHandler*)
 309       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 310     if(aodHandler->GetExtensions()) {
 311
 312       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 313       AliAODEvent *aodFromExt = ext->GetAOD();
 314
 315
 316       switch(fDecChannel){
 317       case 0:
 318         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 319         break;
 320       case 1:
 321         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 322         break;
 323       case 2:
 324         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 325         break;
 326       case 3:
 327         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 328         break;
 329       case 4:
 330         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 331         break;
 332       case 5:
 333         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 334         break;
 335       }
 336     }
 337   } else {
 338     switch(fDecChannel){
 339     case 0:
 340       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 341       break;
 342     case 1:
 343       arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 344       break;
 345     case 2:
 346       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
 347       break;
 348     case 3:
 349       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 350       break;
 351     case 4:
 352       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm4Prong");
 353       break;
 354     case 5:
 355       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 356       break;
 357     }
 358   }
 359   if(!arrayProng) {
 360     printf("AliAnalysisTaskSESignificance::UserExec: branch not found!\n");
 361     return;
 362   }
 363
 364   // fix for temporary bug in ESDfilter
 365   // the AODs with null vertex pointer didn't pass the PhysSel
 366   if(!aod->GetPrimaryVertex()) return;
 367
 368   TClonesArray *arrayMC=0;
 369   AliAODMCHeader *mcHeader=0;
 370
 371   // load MC particles
 372   if(fReadMC){
 373
 374     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 375     if(!arrayMC) {
 376       printf("AliAnalysisTaskSESignificance::UserExec: MC particles branch not found!\n");
 377       //    return;
 378     }
 379
 380     // load MC header
 381     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 382     if(!mcHeader) {
 383       printf("AliAnalysisTaskSESignificance::UserExec: MC header branch not found!\n");
 384       return;
 385     }
 386   }
 387
 388   fHistNEvents->Fill(0); // count event
 389
 390   AliAODRecoDecayHF *d=0;
 391   Int_t nprongs=1;
 392   Int_t *pdgdaughters=0x0;
 393   Int_t prongpdg=211;
 394
 395
 396   switch(fDecChannel){
 397   case 0:
 398     //D+
 399     pdgdaughters =new Int_t[3];
 400     pdgdaughters[0]=211;//pi
 401     pdgdaughters[1]=321;//K
 402     pdgdaughters[2]=211;//pi
 403     nprongs=3;
 404     prongpdg=411;
 405     break;
 406   case 1:
 407     //D0
 408     pdgdaughters =new Int_t[2];
 409     pdgdaughters[0]=211;//pi
 410     pdgdaughters[1]=321;//K
 411     nprongs=2;
 412     prongpdg=421;
 413     break;
 414   case 2:
 415     //D*
 416     pdgdaughters =new Int_t[3];
 417     pdgdaughters[1]=211;//pi
 418     pdgdaughters[0]=321;//K
 419     pdgdaughters[2]=211;//pi (soft?)
 420     nprongs=3;
 421     prongpdg=413;
 422     break;
 423   case 3:
 424     //Ds
 425     pdgdaughters =new Int_t[3];
 426     pdgdaughters[0]=211;//pi
 427     pdgdaughters[1]=321;//K
 428     pdgdaughters[2]=321;//K
 429     nprongs=3;
 430     prongpdg=431;
 431     break;
 432   case 4:
 433     //D0
 434     pdgdaughters =new Int_t[4];
 435     pdgdaughters[0]=321;
 436     pdgdaughters[1]=211;
 437     pdgdaughters[2]=211;
 438     pdgdaughters[3]=211;
 439     nprongs=4;
 440     prongpdg=421;
 441     break;
 442   case 5:
 443     //Lambda_c
 444     pdgdaughters =new Int_t[3];
 445     //check order
 446     pdgdaughters[0]=2212;//p
 447     pdgdaughters[1]=321;//K
 448     pdgdaughters[2]=211;//pi
 449     nprongs=3;
 450     prongpdg=4122;
 451     fReadMC=kFALSE;
 452     printf("Warning! For the moment it is not allowed to use MC for the Lambdac\n");
 453     break;
 454   }
 455   Int_t prongPdgPlus=prongpdg;
 456   Int_t prongPdgMinus=TMath::Abs(prongPdgPlus)*(-1);
 457
 458   Int_t nHistpermv=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
 459   Int_t nProng = arrayProng->GetEntriesFast();
 460   if(fDebug>1) printf("Number of D2H: %d\n",nProng);
 461
 462   if(!fRDCuts->IsEventSelected(aod)){
 463     fHistNEvents->Fill(1);
 464     return;
 465   }
 466
 467   for (Int_t iProng = 0; iProng < nProng; iProng++) {
 468     d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
 469     Double_t invMass=d->InvMass(nprongs,(UInt_t*)pdgdaughters);
 470     Double_t invMassC=0;
 471
 472     Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel);
 473
 474     if(isSelected) {
 475       fHistNEvents->Fill(2); // count selected with loosest cuts
 476       if(fDebug>1) printf("+++++++Is Selected\n");
 477
 478       const Int_t nvars=fRDCuts->GetNVarsForOpt();
 479       Float_t *vars = new Float_t[nvars];
 480       Int_t nVals=0;
 481
 482       fRDCuts->GetCutVarsForOpt(d,vars,nvars,pdgdaughters);
 483       Int_t ptbin=fRDCuts->PtBin(d->Pt());
 484       if(ptbin==-1) continue;
 485       TString mdvname=Form("multiDimVectorPtBin%d",ptbin);
 486       ULong64_t *addresses = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGlobalAddressesAboveCuts(vars,(Float_t)d->Pt(),nVals);
 487       if(fDebug>1)printf("nvals = %d\n",nVals);
 488       for(Int_t ivals=0;ivals<nVals;ivals++){
 489         if(addresses[ivals]>=((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetNTotCells()){
 490           if (fDebug>1) printf("Overflow!!\n");
 491           return;
 492         }
 493         fMassHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
 494         fHistNEvents->Fill(3);
 495         if(fReadMC){
 496           Int_t lab=-1;
 497           lab = d->MatchToMC(prongPdgPlus,arrayMC,nprongs,pdgdaughters);
 498           if(lab>=0){ //signal
 499             AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(lab);
 500             Int_t pdgMC = dMC->GetPdgCode();
 501
 502             if(pdgMC==prongPdgPlus) fSigHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
 503             else fRflHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
 504           }
 505           else{ //background
 506             fBkgHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
 507           }
 508
 509         }
 510       }
 511
 512
 513       if(fDecChannel==AliAnalysisTaskSESignificance::kD0toKpi || fDecChannel==AliAnalysisTaskSESignificance::kLambdactopKpi){ //repeat the cycle to checks cuts passed as D0bar
 514         if(fDecChannel==AliAnalysisTaskSESignificance::kD0toKpi){
 515           pdgdaughters[0]=321;
 516           pdgdaughters[1]=211;
 517         }
 518         if(fDecChannel==AliAnalysisTaskSESignificance::kLambdactopKpi){
 519           pdgdaughters[0]=211;
 520           pdgdaughters[1]=321;
 521           pdgdaughters[2]=2212;
 522         }
 523
 524         invMassC=d->InvMass(nprongs,(UInt_t*)pdgdaughters);
 525         fRDCuts->GetCutVarsForOpt(d,vars,nvars,pdgdaughters);
 526         addresses = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGlobalAddressesAboveCuts(vars,(Float_t)d->Pt(),nVals);
 527         for(Int_t ivals=0;ivals<nVals;ivals++){
 528           fMassHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
 529           if(fReadMC){
 530             Int_t lab=-1;
 531             lab = d->MatchToMC(prongPdgMinus,arrayMC,nprongs,pdgdaughters);
 532             if(lab>=0){ //signal
 533               AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(lab);
 534               Int_t pdgMC = dMC->GetPdgCode();
 535               if(pdgMC==prongPdgMinus) fSigHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
 536               else fRflHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
 537             }
 538             else{ //background
 539               fBkgHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
 540             }
 541
 542           }
 543         }
 544       }
 545
 546     }// end if selected
 547
 548
 549   }
 550
 551   PostData(1,fOutput);
 552   return;
 553 }
 554
 555
 556 //________________________________________________________________________
 557 void AliAnalysisTaskSESignificance::Terminate(Option_t */*option*/)
 558 {
 559   // Terminate analysis
 560   //
 561   if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
 562
 563
 564   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 565   if (!fOutput) {
 566     printf("ERROR: fOutput not available\n");
 567     return;
 568   }
 569
 570   fCutList = dynamic_cast<TList*> (GetOutputData(2));
 571   if (!fCutList) {
 572     printf("ERROR: fCutList not available\n");
 573     return;
 574   }
 575
 576   AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
 577   if (!mdvtmp){
 578     cout<<"multidimvec not found in TList"<<endl;
 579     fCutList->ls();
 580     return;
 581   }
 582   Int_t nHist=mdvtmp->GetNTotCells();
 583   TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
 584   Bool_t drawn=kFALSE;
 585   for(Int_t i=0;i<nHist;i++){
 586
 587     TString hisname;
 588     TString signame;
 589     TString bkgname;
 590     TString rflname;
 591
 592     hisname.Form("hMass_%d",i);
 593     signame.Form("hSig_%d",i);
 594     bkgname.Form("hBkg_%d",i);
 595     rflname.Form("hRfl_%d",i);
 596
 597     fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 598
 599     if (!drawn && fMassHist[i]->GetEntries() > 0){
 600       c1->cd();
 601       fMassHist[i]->Draw();
 602       drawn=kTRUE;
 603     }
 604
 605     if(fReadMC){
 606       fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
 607       fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
 608       fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
 609     }
 610
 611   }
 612
 613
 614
 615
 616   return;
 617 }
 618 //-------------------------------------------
 619