PWG3/vertexingHF/AliAnalysisTaskSEDs.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, 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 //  Analysis task to produce Ds candidates mass spectra          //
  21 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  22 //                                                               //
  23 ///////////////////////////////////////////////////////////////////
  24
  25 #include <TClonesArray.h>
  26 #include <TNtuple.h>
  27 #include <TList.h>
  28 #include <TString.h>
  29 #include <TH1F.h>
  30 #include <TMath.h>
  31 #include <TDatabasePDG.h>
  32
  33 #include "AliAnalysisManager.h"
  34 #include "AliAODHandler.h"
  35 #include "AliAODEvent.h"
  36 #include "AliAODVertex.h"
  37 #include "AliAODTrack.h"
  38 #include "AliAODMCHeader.h"
  39 #include "AliAODMCParticle.h"
  40 #include "AliAODRecoDecayHF3Prong.h"
  41 #include "AliAnalysisVertexingHF.h"
  42 #include "AliRDHFCutsDstoKKpi.h"
  43 #include "AliAnalysisTaskSE.h"
  44 #include "AliAnalysisTaskSEDs.h"
  45 #include "AliNormalizationCounter.h"
  46
  47 ClassImp(AliAnalysisTaskSEDs)
  48
  49
  50 //________________________________________________________________________
  51 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
  52   AliAnalysisTaskSE(),
  53   fOutput(0),
  54   fHistNEvents(0),
  55   fPtVsMass(0),
  56   fYVsPt(0),
  57   fYVsPtSig(0),
  58   fReadMC(kFALSE),
  59   fDoCutVarHistos(kTRUE),
  60   fUseSelectionBit(kFALSE),
  61   fNPtBins(0),
  62   fListCuts(0),
  63   fMassRange(0.4),
  64   fMassBinSize(0.002),
  65   fCounter(0),
  66   fProdCuts(0),
  67   fAnalysisCuts(0)
  68 {
  69   // Default constructor
  70 }
  71
  72 //________________________________________________________________________
  73 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name, AliRDHFCutsDstoKKpi* productioncuts, AliRDHFCutsDstoKKpi* analysiscuts):
  74   AliAnalysisTaskSE(name),
  75   fOutput(0),
  76   fHistNEvents(0),
  77   fPtVsMass(0),
  78   fYVsPt(0),
  79   fYVsPtSig(0),
  80   fReadMC(kFALSE),
  81   fDoCutVarHistos(kTRUE),
  82   fUseSelectionBit(kFALSE),
  83   fNPtBins(0),
  84   fListCuts(0),
  85   fMassRange(0.4),
  86   fMassBinSize(0.002),
  87   fCounter(0),
  88   fProdCuts(productioncuts),
  89   fAnalysisCuts(analysiscuts)
  90 {
  91   // Default constructor
  92   // Output slot #1 writes into a TList container
  93   Int_t nptbins=fAnalysisCuts->GetNPtBins();
  94   Float_t *ptlim=fAnalysisCuts->GetPtBinLimits();
  95   SetPtBins(nptbins,ptlim);
  96
  97   DefineOutput(1,TList::Class());  //My private output
  98
  99   DefineOutput(2,TList::Class());
 100
 101   DefineOutput(3,AliNormalizationCounter::Class());
 102 }
 103
 104 //________________________________________________________________________
 105 void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Float_t* lim){
 106   // define pt bins for analysis
 107   if(n>kMaxPtBins){
 108     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
 109     fNPtBins=kMaxPtBins;
 110     fPtLimits[0]=0.;
 111     fPtLimits[1]=1.;
 112     fPtLimits[2]=3.;
 113     fPtLimits[3]=5.;
 114     fPtLimits[4]=10.;
 115     for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
 116   }else{
 117     fNPtBins=n;
 118     for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
 119     for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
 120   }
 121   if(fDebug > 1){
 122     printf("Number of Pt bins = %d\n",fNPtBins);
 123     for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);
 124   }
 125 }
 126 //________________________________________________________________________
 127 AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()
 128 {
 129   // Destructor
 130   if (fOutput) {
 131     delete fOutput;
 132     fOutput = 0;
 133   }
 134
 135   if(fHistNEvents){
 136     delete fHistNEvents;
 137     fHistNEvents=0;
 138   }
 139
 140   if (fListCuts) {
 141     delete fListCuts;
 142     fListCuts = 0;
 143   }
 144
 145   for(Int_t i=0;i<4*fNPtBins;i++){
 146
 147     if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
 148     if(fMassHistPhi[i]){ delete fMassHistPhi[i]; fMassHistPhi[i]=0;}
 149     if(fMassHistK0st[i]){ delete fMassHistK0st[i]; fMassHistK0st[i]=0;}
 150     if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
 151     if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
 152     if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
 153     if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
 154     if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
 155     if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
 156     if(fPtProng0Hist[i]){ delete fPtProng0Hist[i]; fPtProng0Hist[i]=0;}
 157     if(fPtProng1Hist[i]){ delete fPtProng1Hist[i]; fPtProng1Hist[i]=0;}
 158     if(fPtProng2Hist[i]){ delete fPtProng2Hist[i]; fPtProng2Hist[i]=0;}
 159     if(fDalitz[i]){ delete fDalitz[i]; fDalitz[i]=0;}
 160     if(fDalitzPhi[i]){ delete fDalitzPhi[i]; fDalitzPhi[i]=0;}
 161     if(fDalitzK0st[i]){ delete fDalitzK0st[i]; fDalitzK0st[i]=0;}
 162
 163   }
 164
 165   if(fPtVsMass){
 166     delete fPtVsMass;
 167     fPtVsMass=0;
 168   }
 169   if(fYVsPt){
 170     delete fYVsPt;
 171     fYVsPt=0;
 172   }
 173   if(fYVsPtSig){
 174     delete fYVsPtSig;
 175     fYVsPtSig=0;
 176   }
 177   if(fCounter){
 178     delete fCounter;
 179     fCounter = 0;
 180   }
 181
 182   if (fProdCuts) {
 183     delete fProdCuts;
 184     fProdCuts = 0;
 185   }
 186   if (fAnalysisCuts) {
 187     delete fAnalysisCuts;
 188     fAnalysisCuts = 0;
 189   }
 190 }
 191
 192 //________________________________________________________________________
 193 void AliAnalysisTaskSEDs::Init()
 194 {
 195   // Initialization
 196
 197   if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
 198
 199   fListCuts=new TList();
 200   fListCuts->SetOwner();
 201   fListCuts->SetName("CutObjects");
 202
 203   AliRDHFCutsDstoKKpi *production = new AliRDHFCutsDstoKKpi(*fProdCuts);
 204   production->SetName("ProductionCuts");
 205   AliRDHFCutsDstoKKpi *analysis = new AliRDHFCutsDstoKKpi(*fAnalysisCuts);
 206   analysis->SetName("AnalysisCuts");
 207
 208   fListCuts->Add(production);
 209   fListCuts->Add(analysis);
 210   PostData(2,fListCuts);
 211   return;
 212 }
 213
 214 //________________________________________________________________________
 215 void AliAnalysisTaskSEDs::UserCreateOutputObjects()
 216 {
 217   // Create the output container
 218   //
 219   if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
 220
 221   // Several histograms are more conveniently managed in a TList
 222   fOutput = new TList();
 223   fOutput->SetOwner();
 224   fOutput->SetName("OutputHistos");
 225
 226   Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
 227   Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
 228   if(nInvMassBins%2==1) nInvMassBins++;
 229   // Double_t minMass=massDs-1.0*nInvMassBins*fMassBinSize;
 230   Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
 231   //  Double_t maxMass=massDs+1.0*nInvMassBins*fMassBinSize;
 232   Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
 233
 234   TString hisname;
 235   TString htype;
 236   Int_t index;
 237   for(Int_t iType=0; iType<4; iType++){
 238     for(Int_t i=0;i<fNPtBins;i++){
 239       if(iType==0){
 240         htype="All";
 241         index=GetHistoIndex(i);
 242       }else if(iType==1){
 243         htype="Sig";
 244         index=GetSignalHistoIndex(i);
 245       }else if(iType==2){
 246         htype="Bkg";
 247         index=GetBackgroundHistoIndex(i);
 248       }else{
 249         htype="ReflSig";
 250         index=GetReflSignalHistoIndex(i);
 251       }
 252       hisname.Form("hMass%sPt%d",htype.Data(),i);
 253       fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 254       fMassHist[index]->Sumw2();
 255       hisname.Form("hMass%sPt%dphi",htype.Data(),i);
 256       fMassHistPhi[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 257       fMassHistPhi[index]->Sumw2();
 258       hisname.Form("hMass%sPt%dk0st",htype.Data(),i);
 259       fMassHistK0st[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 260       fMassHistK0st[index]->Sumw2();
 261       hisname.Form("hCosP%sPt%d",htype.Data(),i);
 262       fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 263       fCosPHist[index]->Sumw2();
 264       hisname.Form("hDLen%sPt%d",htype.Data(),i);
 265       fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 266       fDLenHist[index]->Sumw2();
 267       hisname.Form("hSumd02%sPt%d",htype.Data(),i);
 268       fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 269       fSumd02Hist[index]->Sumw2();
 270       hisname.Form("hSigVert%sPt%d",htype.Data(),i);
 271       fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 272       fSigVertHist[index]->Sumw2();
 273       hisname.Form("hPtMax%sPt%d",htype.Data(),i);
 274       fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
 275       fPtMaxHist[index]->Sumw2();
 276       hisname.Form("hPtCand%sPt%d",htype.Data(),i);
 277       fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
 278       fPtCandHist[index]->Sumw2();
 279       hisname.Form("hDCA%sPt%d",htype.Data(),i);
 280       fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 281       fDCAHist[index]->Sumw2();
 282       hisname.Form("hPtProng0%sPt%d",htype.Data(),i);
 283       fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 284       fPtProng0Hist[index]->Sumw2();
 285       hisname.Form("hPtProng1%sPt%d",htype.Data(),i);
 286       fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 287       fPtProng1Hist[index]->Sumw2();
 288       hisname.Form("hPtProng2%sPt%d",htype.Data(),i);
 289       fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 290       fPtProng2Hist[index]->Sumw2();
 291       hisname.Form("hDalitz%sPt%d",htype.Data(),i);
 292       fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 293       fDalitz[index]->Sumw2();
 294       hisname.Form("hDalitz%sPt%dphi",htype.Data(),i);
 295       fDalitzPhi[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 296       fDalitzPhi[index]->Sumw2();
 297       hisname.Form("hDalitz%sPt%dk0st",htype.Data(),i);
 298       fDalitzK0st[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 299       fDalitzK0st[index]->Sumw2();
 300     }
 301   }
 302
 303   for(Int_t i=0; i<4*fNPtBins; i++){
 304     fOutput->Add(fMassHist[i]);
 305     fOutput->Add(fMassHistPhi[i]);
 306     fOutput->Add(fMassHistK0st[i]);
 307     fOutput->Add(fCosPHist[i]);
 308     fOutput->Add(fDLenHist[i]);
 309     fOutput->Add(fSumd02Hist[i]);
 310     fOutput->Add(fSigVertHist[i]);
 311     fOutput->Add(fPtMaxHist[i]);
 312     fOutput->Add(fPtCandHist[i]);
 313     fOutput->Add(fDCAHist[i]);
 314     fOutput->Add(fPtProng0Hist[i]);
 315     fOutput->Add(fPtProng1Hist[i]);
 316     fOutput->Add(fPtProng2Hist[i]);
 317     fOutput->Add(fDalitz[i]);
 318     fOutput->Add(fDalitzPhi[i]);
 319     fOutput->Add(fDalitzK0st[i]);
 320   }
 321
 322   fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",64,-0.5,63.5);
 323   fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",64,-0.5,63.5);
 324   fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",64,-0.5,63.5);
 325   fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",64,-0.5,63.5);
 326   for(Int_t i=0;i<4;i++){
 327     fChanHist[i]->Sumw2();
 328     fChanHist[i]->SetMinimum(0);
 329     fOutput->Add(fChanHist[i]);
 330   }
 331
 332   fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
 333   fHistNEvents->Sumw2();
 334   fHistNEvents->SetMinimum(0);
 335   fOutput->Add(fHistNEvents);
 336
 337   fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nInvMassBins,minMass,maxMass,40,0.,20.);
 338   fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
 339   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
 340
 341   fOutput->Add(fPtVsMass);
 342   fOutput->Add(fYVsPt);
 343   fOutput->Add(fYVsPtSig);
 344
 345   //Counter for Normalization
 346   fCounter = new AliNormalizationCounter("NormalizationCounter");
 347
 348   PostData(1,fOutput);
 349   PostData(3,fCounter);
 350
 351   return;
 352 }
 353
 354 //________________________________________________________________________
 355 void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
 356 {
 357   // Ds selection for current event, fill mass histos and selecetion variable histo
 358   // separate signal and backgound if fReadMC is activated
 359
 360   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 361
 362
 363   TClonesArray *array3Prong = 0;
 364   if(!aod && AODEvent() && IsStandardAOD()) {
 365     // In case there is an AOD handler writing a standard AOD, use the AOD
 366     // event in memory rather than the input (ESD) event.
 367     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 368     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 369     // have to taken from the AOD event hold by the AliAODExtension
 370     AliAODHandler* aodHandler = (AliAODHandler*)
 371       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 372     if(aodHandler->GetExtensions()) {
 373       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 374       AliAODEvent *aodFromExt = ext->GetAOD();
 375       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 376     }
 377   } else if(aod) {
 378     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 379   }
 380
 381   if(!aod || !array3Prong) {
 382     printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
 383     return;
 384   }
 385
 386
 387   // fix for temporary bug in ESDfilter
 388   // the AODs with null vertex pointer didn't pass the PhysSel
 389   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 390
 391
 392   fHistNEvents->Fill(0); // count event
 393   // Post the data already here
 394   PostData(1,fOutput);
 395
 396   fCounter->StoreEvent(aod,fReadMC);
 397
 398   TClonesArray *arrayMC=0;
 399   AliAODMCHeader *mcHeader=0;
 400
 401   // AOD primary vertex
 402   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 403   //    vtx1->Print();
 404
 405   // load MC particles
 406   if(fReadMC){
 407
 408     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 409     if(!arrayMC) {
 410       printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
 411       return;
 412     }
 413
 414   // load MC header
 415     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 416     if(!mcHeader) {
 417       printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
 418       return;
 419     }
 420   }
 421
 422   Int_t n3Prong = array3Prong->GetEntriesFast();
 423   if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
 424
 425
 426   Int_t pdgDstoKKpi[3]={321,321,211};
 427   Int_t nSelectedloose=0;
 428   Int_t nSelectedtight=0;
 429
 430   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 431     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 432
 433     if(fUseSelectionBit && !(d->HasSelectionBit(AliRDHFCuts::kDsCuts))) continue;
 434
 435     Bool_t unsetvtx=kFALSE;
 436     if(!d->GetOwnPrimaryVtx()){
 437       d->SetOwnPrimaryVtx(vtx1);
 438       unsetvtx=kTRUE;
 439     }
 440
 441     Double_t ptCand = d->Pt();
 442     Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
 443     Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate);
 444     Double_t rapid=d->YDs();
 445     fYVsPt->Fill(ptCand,rapid);
 446
 447     Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
 448     if(isFidAcc){
 449       nSelectedloose++;
 450       if(retCodeAnalysisCuts>0)nSelectedtight++;
 451     }
 452     if(retCodeAnalysisCuts<=0) continue;
 453     if(!isFidAcc) continue;
 454
 455     Int_t index=GetHistoIndex(iPtBin);
 456     fPtCandHist[index]->Fill(ptCand);
 457
 458     Int_t isKKpi=retCodeAnalysisCuts&1;
 459     Int_t ispiKK=retCodeAnalysisCuts&2;
 460     Int_t isPhiKKpi=retCodeAnalysisCuts&4;
 461     Int_t isPhipiKK=retCodeAnalysisCuts&8;
 462     Int_t isK0starKKpi=retCodeAnalysisCuts&16;
 463     Int_t isK0starpiKK=retCodeAnalysisCuts&32;
 464
 465     fChanHist[0]->Fill(retCodeAnalysisCuts);
 466
 467     Int_t indexMCKKpi=-1;
 468     Int_t indexMCpiKK=-1;
 469     Int_t labDs=-1;
 470     if(fReadMC){
 471       labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
 472       if(labDs>=0){
 473         Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 474         AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
 475         Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 476
 477         if(isKKpi){
 478           if(pdgCode0==321) {
 479             indexMCKKpi=GetSignalHistoIndex(iPtBin);
 480             fYVsPtSig->Fill(ptCand,rapid);
 481             fChanHist[1]->Fill(retCodeAnalysisCuts);
 482           }else{
 483             indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
 484             fChanHist[3]->Fill(retCodeAnalysisCuts);
 485           }
 486         }
 487         if(ispiKK){
 488           if(pdgCode0==211) {
 489             indexMCpiKK=GetSignalHistoIndex(iPtBin);
 490             fYVsPtSig->Fill(ptCand,rapid);
 491             fChanHist[1]->Fill(retCodeAnalysisCuts);
 492           }else{
 493             indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
 494             fChanHist[3]->Fill(retCodeAnalysisCuts);
 495           }
 496         }
 497       }else{
 498         indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
 499         indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
 500         fChanHist[2]->Fill(retCodeAnalysisCuts);
 501       }
 502     }
 503
 504     if(isKKpi){
 505       Double_t invMass=d->InvMassDsKKpi();
 506       fMassHist[index]->Fill(invMass);
 507       fPtVsMass->Fill(invMass,ptCand);
 508       if(isPhiKKpi) fMassHistPhi[index]->Fill(invMass);
 509       if(isK0starKKpi) fMassHistK0st[index]->Fill(invMass);
 510       if(fReadMC  && indexMCKKpi!=-1){
 511         fMassHist[indexMCKKpi]->Fill(invMass);
 512         if(isPhiKKpi) fMassHistPhi[indexMCKKpi]->Fill(invMass);
 513         if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass);
 514       }
 515     }
 516     if(ispiKK){
 517       Double_t invMass=d->InvMassDspiKK();
 518       fMassHist[index]->Fill(invMass);
 519       fPtVsMass->Fill(invMass,ptCand);
 520       if(isPhipiKK) fMassHistPhi[index]->Fill(invMass);
 521       if(isK0starpiKK) fMassHistK0st[index]->Fill(invMass);
 522       if(fReadMC  && indexMCpiKK!=-1){
 523         fMassHist[indexMCpiKK]->Fill(invMass);
 524         if(isPhipiKK) fMassHistPhi[indexMCpiKK]->Fill(invMass);
 525         if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass);
 526       }
 527     }
 528
 529     if(fDoCutVarHistos){
 530       Double_t dlen=d->DecayLength();
 531       Double_t cosp=d->CosPointingAngle();
 532       Double_t pt0=d->PtProng(0);
 533       Double_t pt1=d->PtProng(1);
 534       Double_t pt2=d->PtProng(2);
 535       Double_t sigvert=d->GetSigmaVert();
 536       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 537       Double_t dca=d->GetDCA();
 538       Double_t ptmax=0;
 539       for(Int_t i=0;i<3;i++){
 540         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 541       }
 542       fCosPHist[index]->Fill(cosp);
 543       fDLenHist[index]->Fill(dlen);
 544       fSigVertHist[index]->Fill(sigvert);
 545       fSumd02Hist[index]->Fill(sumD02);
 546       fPtMaxHist[index]->Fill(ptmax);
 547       fDCAHist[index]->Fill(dca);
 548       fPtProng0Hist[index]->Fill(pt0);
 549       fPtProng1Hist[index]->Fill(pt1);
 550       fPtProng2Hist[index]->Fill(pt2);
 551       if(isKKpi){
 552         Double_t massKK=d->InvMass2Prongs(0,1,321,321);
 553         Double_t massKp=d->InvMass2Prongs(1,2,321,211);
 554         fDalitz[index]->Fill(massKK,massKp);
 555         if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
 556         if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
 557         if(fReadMC && indexMCKKpi!=-1){
 558           fDalitz[indexMCKKpi]->Fill(massKK,massKp);
 559           if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
 560           if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
 561           fCosPHist[indexMCKKpi]->Fill(cosp);
 562           fDLenHist[indexMCKKpi]->Fill(dlen);
 563           fSigVertHist[indexMCKKpi]->Fill(sigvert);
 564           fSumd02Hist[indexMCKKpi]->Fill(sumD02);
 565           fPtMaxHist[indexMCKKpi]->Fill(ptmax);
 566           fPtCandHist[indexMCKKpi]->Fill(ptCand);
 567           fDCAHist[indexMCKKpi]->Fill(dca);
 568           fPtProng0Hist[indexMCKKpi]->Fill(pt0);
 569           fPtProng1Hist[indexMCKKpi]->Fill(pt1);
 570           fPtProng2Hist[indexMCKKpi]->Fill(pt2);
 571         }
 572       }
 573       if(ispiKK){
 574         Double_t massKK=d->InvMass2Prongs(1,2,321,321);
 575         Double_t massKp=d->InvMass2Prongs(0,1,211,321);
 576         fDalitz[index]->Fill(massKK,massKp);
 577         if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
 578         if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
 579         if(fReadMC && indexMCpiKK!=-1){
 580           fDalitz[indexMCpiKK]->Fill(massKK,massKp);
 581           if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
 582           if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
 583           fCosPHist[indexMCpiKK]->Fill(cosp);
 584           fDLenHist[indexMCpiKK]->Fill(dlen);
 585           fSigVertHist[indexMCpiKK]->Fill(sigvert);
 586           fSumd02Hist[indexMCpiKK]->Fill(sumD02);
 587           fPtMaxHist[indexMCpiKK]->Fill(ptmax);
 588           fPtCandHist[indexMCpiKK]->Fill(ptCand);
 589           fDCAHist[indexMCpiKK]->Fill(dca);
 590           fPtProng0Hist[indexMCpiKK]->Fill(pt0);
 591           fPtProng1Hist[indexMCpiKK]->Fill(pt1);
 592           fPtProng2Hist[indexMCpiKK]->Fill(pt2);
 593         }
 594       }
 595
 596     }
 597
 598     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 599   }
 600
 601   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
 602   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
 603
 604   PostData(1,fOutput);
 605   PostData(3,fCounter);
 606
 607   return;
 608 }
 609
 610 //_________________________________________________________________
 611
 612 void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
 613 {
 614   // Terminate analysis
 615   //
 616   if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
 617   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 618   if (!fOutput) {
 619     printf("ERROR: fOutput not available\n");
 620     return;
 621   }
 622   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
 623   if(fHistNEvents){
 624     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
 625   }else{
 626     printf("ERROR: fHistNEvents not available\n");
 627     return;
 628   }
 629   return;
 630 }
 631