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