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   fHistNEvents->Fill(0); // count event
 305   // Post the data already here
 306   PostData(1,fOutput);
 307
 308
 309   TClonesArray *array3Prong = 0;
 310   if(!aod && AODEvent() && IsStandardAOD()) {
 311     // In case there is an AOD handler writing a standard AOD, use the AOD
 312     // event in memory rather than the input (ESD) event.
 313     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 314     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 315     // have to taken from the AOD event hold by the AliAODExtension
 316     AliAODHandler* aodHandler = (AliAODHandler*)
 317       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 318     if(aodHandler->GetExtensions()) {
 319       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 320       AliAODEvent *aodFromExt = ext->GetAOD();
 321       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 322     }
 323   } else {
 324     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 325   }
 326
 327   if(!array3Prong) {
 328     printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
 329     return;
 330   }
 331
 332
 333   TClonesArray *arrayMC=0;
 334   AliAODMCHeader *mcHeader=0;
 335
 336   // AOD primary vertex
 337   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 338   //    vtx1->Print();
 339
 340   // load MC particles
 341   if(fReadMC){
 342
 343     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 344     if(!arrayMC) {
 345       printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
 346       return;
 347     }
 348
 349   // load MC header
 350     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 351     if(!mcHeader) {
 352       printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
 353       return;
 354     }
 355   }
 356
 357   Int_t n3Prong = array3Prong->GetEntriesFast();
 358   if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
 359
 360
 361   Int_t pdgDstoKKpi[3]={321,321,211};
 362   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 363     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 364
 365
 366     Bool_t unsetvtx=kFALSE;
 367     if(!d->GetOwnPrimaryVtx()){
 368       d->SetOwnPrimaryVtx(vtx1);
 369       unsetvtx=kTRUE;
 370     }
 371
 372     Int_t retCodeProductionCuts=fProdCuts->IsSelected(d,AliRDHFCuts::kCandidate);
 373     if(retCodeProductionCuts>0){
 374       Int_t isKKpi=retCodeProductionCuts&1;
 375       Int_t ispiKK=retCodeProductionCuts&2;
 376       Int_t isPhi=retCodeProductionCuts&4;
 377       Int_t isK0star=retCodeProductionCuts&8;
 378       Double_t ptCand = d->Pt();
 379       Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
 380       Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate);
 381       Int_t isKKpiAC=retCodeAnalysisCuts&1;
 382       Int_t ispiKKAC=retCodeAnalysisCuts&2;
 383       Int_t isPhiAC=retCodeAnalysisCuts&4;
 384       Int_t isK0starAC=retCodeAnalysisCuts&8;
 385
 386       Int_t labDs=-1;
 387       if(fReadMC){
 388         labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
 389       }
 390
 391       Double_t dlen=d->DecayLength();
 392       Double_t cosp=d->CosPointingAngle();
 393       Double_t pt0=d->PtProng(0);
 394       Double_t pt1=d->PtProng(1);
 395       Double_t pt2=d->PtProng(2);
 396       Double_t sigvert=d->GetSigmaVert();
 397       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 398       Double_t dca=d->GetDCA();
 399       Double_t ptmax=0;
 400       for(Int_t i=0;i<3;i++){
 401         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 402       }
 403
 404       Int_t index=GetHistoIndex(iPtBin);
 405       Int_t type=0;
 406       if(isKKpi) type+=1;
 407       if(ispiKK) type+=2;
 408       Int_t typeAC=0;
 409       if(isKKpiAC) typeAC+=1;
 410       if(ispiKKAC) typeAC+=2;
 411       fCosPHist[index]->Fill(cosp);
 412       fDLenHist[index]->Fill(dlen);
 413       fSigVertHist[index]->Fill(sigvert);
 414       fSumd02Hist[index]->Fill(sumD02);
 415       fPtMaxHist[index]->Fill(ptmax);
 416       fPtCandHist[index]->Fill(ptCand);
 417       fDCAHist[index]->Fill(dca);
 418       fChanHist[0]->Fill(type);
 419       fPtProng0Hist[index]->Fill(pt0);
 420       fPtProng1Hist[index]->Fill(pt1);
 421       fPtProng2Hist[index]->Fill(pt2);
 422
 423       if(retCodeAnalysisCuts>0) fChanHistCuts[0]->Fill(typeAC);
 424       if(fReadMC){
 425         if(labDs>=0) {
 426           index=GetSignalHistoIndex(iPtBin);
 427           fChanHist[1]->Fill(type);
 428           if(retCodeAnalysisCuts>0) fChanHistCuts[1]->Fill(typeAC);
 429         }else{
 430           index=GetBackgroundHistoIndex(iPtBin);
 431           fChanHist[2]->Fill(type);
 432           if(retCodeAnalysisCuts>0) fChanHistCuts[2]->Fill(typeAC);
 433         }
 434       }
 435       if(isKKpi){
 436         index=GetHistoIndex(iPtBin);
 437         Double_t invMass=d->InvMassDsKKpi();
 438         fMassHist[index]->Fill(invMass);
 439         if(isPhi) fMassHistPhi[index]->Fill(invMass);
 440         if(isK0star) fMassHistK0st[index]->Fill(invMass);
 441         Double_t massKK=d->InvMass2Prongs(0,1,321,321);
 442         Double_t massKp=d->InvMass2Prongs(1,2,321,211);
 443         fDalitz[index]->Fill(massKK,massKp);
 444         if(isPhi) fDalitzPhi[index]->Fill(massKK,massKp);
 445         if(isK0star) fDalitzK0st[index]->Fill(massKK,massKp);
 446         if(retCodeAnalysisCuts>0 && isKKpiAC){
 447           fMassHistCuts[index]->Fill(invMass);
 448           if(isPhiAC) fMassHistCutsPhi[index]->Fill(invMass);
 449           if(isK0starAC) fMassHistCutsK0st[index]->Fill(invMass);
 450         }
 451         if(fReadMC){
 452           Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 453           AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 454           Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 455           if(labDs>=0){
 456             if(pdgCode0==321) {
 457               index=GetSignalHistoIndex(iPtBin);
 458             }else{
 459               index=GetReflSignalHistoIndex(iPtBin);
 460             }
 461           }else{
 462             index=GetBackgroundHistoIndex(iPtBin);
 463           }
 464           fMassHist[index]->Fill(invMass);
 465           if(isPhi) fMassHistPhi[index]->Fill(invMass);
 466           if(isK0star) fMassHistK0st[index]->Fill(invMass);
 467           fCosPHist[index]->Fill(cosp);
 468           fDLenHist[index]->Fill(dlen);
 469           fSigVertHist[index]->Fill(sigvert);
 470           fSumd02Hist[index]->Fill(sumD02);
 471           fPtMaxHist[index]->Fill(ptmax);
 472           fPtCandHist[index]->Fill(ptCand);
 473           fDCAHist[index]->Fill(dca);
 474           fPtProng0Hist[index]->Fill(pt0);
 475           fPtProng1Hist[index]->Fill(pt1);
 476           fPtProng2Hist[index]->Fill(pt2);
 477           fDalitz[index]->Fill(massKK,massKp);
 478           if(isPhi) fDalitzPhi[index]->Fill(massKK,massKp);
 479           if(isK0star) fDalitzK0st[index]->Fill(massKK,massKp);
 480           if(retCodeAnalysisCuts>0 && isKKpiAC){
 481             fMassHistCuts[index]->Fill(invMass);
 482             if(isPhiAC) fMassHistCutsPhi[index]->Fill(invMass);
 483             if(isK0starAC) fMassHistCutsK0st[index]->Fill(invMass);
 484           }
 485         }
 486       }
 487       if(ispiKK){
 488         index=GetHistoIndex(iPtBin);
 489         Double_t invMass=d->InvMassDspiKK();
 490         fMassHist[index]->Fill(invMass);
 491         if(isPhi) fMassHistPhi[index]->Fill(invMass);
 492         if(isK0star) fMassHistK0st[index]->Fill(invMass);
 493         Double_t massKK=d->InvMass2Prongs(1,2,321,321);
 494         Double_t massKp=d->InvMass2Prongs(0,1,211,321);
 495         fDalitz[index]->Fill(massKK,massKp);
 496         if(isPhi) fDalitzPhi[index]->Fill(massKK,massKp);
 497         if(isK0star) fDalitzK0st[index]->Fill(massKK,massKp);
 498         if(retCodeAnalysisCuts>0 && ispiKKAC){
 499           fMassHistCuts[index]->Fill(invMass);
 500           if(isPhiAC) fMassHistCutsPhi[index]->Fill(invMass);
 501           if(isK0starAC) fMassHistCutsK0st[index]->Fill(invMass);
 502         }
 503         if(fReadMC){
 504           Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 505           AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 506           Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 507           if(labDs>=0) {
 508             if(pdgCode0==211) {
 509               index=GetSignalHistoIndex(iPtBin);
 510             }else{
 511               index=GetReflSignalHistoIndex(iPtBin);
 512             }
 513           }else{
 514             index=GetBackgroundHistoIndex(iPtBin);
 515           }
 516           fMassHist[index]->Fill(invMass);
 517           if(isPhi) fMassHistPhi[index]->Fill(invMass);
 518           if(isK0star) fMassHistK0st[index]->Fill(invMass);
 519           fCosPHist[index]->Fill(cosp);
 520           fDLenHist[index]->Fill(dlen);
 521           fSigVertHist[index]->Fill(sigvert);
 522           fSumd02Hist[index]->Fill(sumD02);
 523           fPtMaxHist[index]->Fill(ptmax);
 524           fPtCandHist[index]->Fill(ptCand);
 525           fDCAHist[index]->Fill(dca);
 526           fPtProng0Hist[index]->Fill(pt0);
 527           fPtProng1Hist[index]->Fill(pt1);
 528           fPtProng2Hist[index]->Fill(pt2);
 529           fDalitz[index]->Fill(massKK,massKp);
 530           if(isPhi) fDalitzPhi[index]->Fill(massKK,massKp);
 531           if(isK0star) fDalitzK0st[index]->Fill(massKK,massKp);
 532           if(retCodeAnalysisCuts>0 && ispiKKAC){
 533             fMassHistCuts[index]->Fill(invMass);
 534             if(isPhiAC) fMassHistCutsPhi[index]->Fill(invMass);
 535             if(isK0starAC) fMassHistCutsK0st[index]->Fill(invMass);
 536           }
 537         }
 538       }
 539     }
 540     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 541   }
 542
 543
 544   PostData(1,fOutput);
 545   return;
 546 }
 547
 548 //_________________________________________________________________
 549
 550 void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
 551 {
 552   // Terminate analysis
 553   //
 554   if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
 555   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 556   if (!fOutput) {
 557     printf("ERROR: fOutput not available\n");
 558     return;
 559   }
 560   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
 561   fChanHist[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAll"));
 562   fChanHist[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSig"));
 563   fChanHist[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkg"));
 564   fChanHist[3] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanReflSig"));
 565   fChanHistCuts[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAllCuts"));
 566   fChanHistCuts[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSigCuts"));
 567   fChanHistCuts[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkgCuts"));
 568   fChanHistCuts[3] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanReflSigCuts"));
 569
 570
 571   TString hisname;
 572   TString htype;
 573   Int_t index;
 574   for(Int_t iType=0; iType<4; iType++){
 575     for(Int_t i=0;i<fNPtBins;i++){
 576       if(iType==0){
 577         htype="All";
 578         index=GetHistoIndex(i);
 579       }else if(iType==1){
 580         htype="Sig";
 581         index=GetSignalHistoIndex(i);
 582       }else if(iType==2){
 583         htype="Bkg";
 584         index=GetBackgroundHistoIndex(i);
 585       }else{
 586         htype="ReflSig";
 587         index=GetReflSignalHistoIndex(i);
 588       }
 589       hisname.Form("hMass%sPt%d",htype.Data(),i);
 590       fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 591       hisname.Form("hMass%sPt%dCuts",htype.Data(),i);
 592       fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 593       hisname.Form("hMass%sPt%dphi",htype.Data(),i);
 594       fMassHistPhi[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 595       hisname.Form("hMass%sPt%dphiCuts",htype.Data(),i);
 596       fMassHistCutsPhi[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 597       hisname.Form("hMass%sPt%dk0st",htype.Data(),i);
 598       fMassHistK0st[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 599       hisname.Form("hMass%sPt%dk0stCuts",htype.Data(),i);
 600       fMassHistCutsK0st[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 601       hisname.Form("hCosP%sPt%d",htype.Data(),i);
 602       fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 603       hisname.Form("hDLen%sPt%d",htype.Data(),i);
 604       fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 605       hisname.Form("hSumd02%sPt%d",htype.Data(),i);
 606       fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 607       hisname.Form("hSigVert%sPt%d",htype.Data(),i);
 608       fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 609       hisname.Form("hPtMax%sPt%d",htype.Data(),i);
 610       fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 611       hisname.Form("hPtCand%sPt%d",htype.Data(),i);
 612       fPtCandHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 613       hisname.Form("hDCA%sPt%d",htype.Data(),i);
 614       fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 615       hisname.Form("hPtProng0%sPt%d",htype.Data(),i);
 616       fPtProng0Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 617       hisname.Form("hPtProng1%sPt%d",htype.Data(),i);
 618       fPtProng1Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 619       hisname.Form("hPtProng2%sPt%d",htype.Data(),i);
 620       fPtProng2Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 621       hisname.Form("hDalitz%sPt%d",htype.Data(),i);
 622       fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
 623       hisname.Form("hDalitz%sPt%dphi",htype.Data(),i);
 624       fDalitzPhi[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
 625       hisname.Form("hDalitz%sPt%dk0st",htype.Data(),i);
 626       fDalitzK0st[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
 627     }
 628   }
 629   return;
 630 }