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