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   Int_t index;
 173   for(Int_t i=0;i<fNPtBins;i++){
 174     index=GetHistoIndex(i);
 175     hisname.Form("hMassAllPt%d",i);
 176     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 177     fMassHist[index]->Sumw2();
 178     hisname.Form("hMassAllPt%dCuts",i);
 179     fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 180     fMassHistCuts[index]->Sumw2();
 181     hisname.Form("hCosPAllPt%d",i);
 182     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 183     fCosPHist[index]->Sumw2();
 184     hisname.Form("hDLenAllPt%d",i);
 185     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 186     fDLenHist[index]->Sumw2();
 187     hisname.Form("hSumd02AllPt%d",i);
 188     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 189     fSumd02Hist[index]->Sumw2();
 190      hisname.Form("hSigVertAllPt%d",i);
 191     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 192     fSigVertHist[index]->Sumw2();
 193     hisname.Form("hPtMaxAllPt%d",i);
 194     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
 195     fPtMaxHist[index]->Sumw2();
 196     hisname.Form("hPtCandAllPt%d",i);
 197     fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
 198     fPtCandHist[index]->Sumw2();
 199      hisname.Form("hDCAAllPt%d",i);
 200     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 201     fDCAHist[index]->Sumw2();
 202     hisname.Form("hPtProng0AllPt%d",i);
 203     fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 204     fPtProng0Hist[index]->Sumw2();
 205     hisname.Form("hPtProng1AllPt%d",i);
 206     fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 207     fPtProng1Hist[index]->Sumw2();
 208     hisname.Form("hPtProng2AllPt%d",i);
 209     fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 210     fPtProng2Hist[index]->Sumw2();
 211
 212     hisname.Form("hDalitzAllPt%d",i);
 213     fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 214     fDalitz[index]->Sumw2();
 215
 216     index=GetSignalHistoIndex(i);
 217     hisname.Form("hMassSigPt%d",i);
 218     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 219     fMassHist[index]->Sumw2();
 220     hisname.Form("hMassSigPt%dCuts",i);
 221     fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 222     fMassHistCuts[index]->Sumw2();
 223     hisname.Form("hCosPSigPt%d",i);
 224     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 225     fCosPHist[index]->Sumw2();
 226     hisname.Form("hDLenSigPt%d",i);
 227     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 228     fDLenHist[index]->Sumw2();
 229     hisname.Form("hSumd02SigPt%d",i);
 230     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 231     fSumd02Hist[index]->Sumw2();
 232     hisname.Form("hSigVertSigPt%d",i);
 233     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 234     fSigVertHist[index]->Sumw2();
 235     hisname.Form("hPtMaxSigPt%d",i);
 236     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 237     fPtMaxHist[index]->Sumw2();
 238     hisname.Form("hPtCandSigPt%d",i);
 239     fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 240     fPtCandHist[index]->Sumw2();
 241     hisname.Form("hDCASigPt%d",i);
 242     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 243     fDCAHist[index]->Sumw2();
 244     hisname.Form("hPtProng0SigPt%d",i);
 245     fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 246     fPtProng0Hist[index]->Sumw2();
 247     hisname.Form("hPtProng1SigPt%d",i);
 248     fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 249     fPtProng1Hist[index]->Sumw2();
 250     hisname.Form("hPtProng2SigPt%d",i);
 251     fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 252     fPtProng2Hist[index]->Sumw2();
 253     hisname.Form("hDalitzSigPt%d",i);
 254     fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 255     fDalitz[index]->Sumw2();
 256
 257     index=GetBackgroundHistoIndex(i);
 258     hisname.Form("hMassBkgPt%d",i);
 259     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 260     fMassHist[index]->Sumw2();
 261     hisname.Form("hMassBkgPt%dCuts",i);
 262     fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 263     fMassHistCuts[index]->Sumw2();
 264     hisname.Form("hCosPBkgPt%d",i);
 265     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 266     fCosPHist[index]->Sumw2();
 267     hisname.Form("hDLenBkgPt%d",i);
 268     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 269     fDLenHist[index]->Sumw2();
 270     hisname.Form("hSumd02BkgPt%d",i);
 271     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 272     fSumd02Hist[index]->Sumw2();
 273     hisname.Form("hSigVertBkgPt%d",i);
 274     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 275     fSigVertHist[index]->Sumw2();
 276     hisname.Form("hPtMaxSigBkg%d",i);
 277     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 278     fPtMaxHist[index]->Sumw2();
 279     hisname.Form("hPtCandBkgPt%d",i);
 280     fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 281     fPtCandHist[index]->Sumw2();
 282     hisname.Form("hDCABkgPt%d",i);
 283     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 284     fDCAHist[index]->Sumw2();
 285     hisname.Form("hPtProng0BkgPt%d",i);
 286     fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 287     fPtProng0Hist[index]->Sumw2();
 288     hisname.Form("hPtProng1BkgPt%d",i);
 289     fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 290     fPtProng1Hist[index]->Sumw2();
 291     hisname.Form("hPtProng2BkgPt%d",i);
 292     fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 293     fPtProng2Hist[index]->Sumw2();
 294     hisname.Form("hDalitzBkgPt%d",i);
 295     fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 296     fDalitz[index]->Sumw2();
 297
 298     index=GetReflSignalHistoIndex(i);
 299     hisname.Form("hMassReflSigPt%d",i);
 300     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
 301     fMassHist[index]->Sumw2();
 302     hisname.Form("hMassReflSigPt%dCuts",i);
 303     fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
 304     fMassHistCuts[index]->Sumw2();
 305     hisname.Form("hCosPReflSigPt%d",i);
 306     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 307     fCosPHist[index]->Sumw2();
 308     hisname.Form("hDLenReflSigPt%d",i);
 309     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 310     fDLenHist[index]->Sumw2();
 311     hisname.Form("hSumd02ReflSigPt%d",i);
 312     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 313     fSumd02Hist[index]->Sumw2();
 314     hisname.Form("hSigVertReflSigPt%d",i);
 315     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 316     fSigVertHist[index]->Sumw2();
 317     hisname.Form("hPtMaxReflSigPt%d",i);
 318     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 319     fPtMaxHist[index]->Sumw2();
 320     hisname.Form("hPtCandReflSigPt%d",i);
 321     fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 322     fPtCandHist[index]->Sumw2();
 323     hisname.Form("hDCAReflSigPt%d",i);
 324     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 325     fDCAHist[index]->Sumw2();
 326     hisname.Form("hPtProng0ReflSigPt%d",i);
 327     fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 328     fPtProng0Hist[index]->Sumw2();
 329     hisname.Form("hPtProng1ReflSigPt%d",i);
 330     fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 331     fPtProng1Hist[index]->Sumw2();
 332     hisname.Form("hPtProng2ReflSigPt%d",i);
 333     fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 334     fPtProng2Hist[index]->Sumw2();
 335     hisname.Form("hDalitzReflSigPt%d",i);
 336     fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 337     fDalitz[index]->Sumw2();
 338   }
 339
 340   for(Int_t i=0; i<4*fNPtBins; i++){
 341     fOutput->Add(fMassHist[i]);
 342     fOutput->Add(fMassHistCuts[i]);
 343     fOutput->Add(fCosPHist[i]);
 344     fOutput->Add(fDLenHist[i]);
 345     fOutput->Add(fSumd02Hist[i]);
 346     fOutput->Add(fSigVertHist[i]);
 347     fOutput->Add(fPtMaxHist[i]);
 348     fOutput->Add(fPtCandHist[i]);
 349     fOutput->Add(fDCAHist[i]);
 350     fOutput->Add(fPtProng0Hist[i]);
 351     fOutput->Add(fPtProng1Hist[i]);
 352     fOutput->Add(fPtProng2Hist[i]);
 353     fOutput->Add(fDalitz[i]);
 354   }
 355
 356   fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",4,-0.5,3.5);
 357   fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",4,-0.5,3.5);
 358   fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",4,-0.5,3.5);
 359   fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",4,-0.5,3.5);
 360   fChanHistCuts[0] = new TH1F("hChanAllCuts", "KKpi and piKK candidates",4,-0.5,3.5);
 361   fChanHistCuts[1] = new TH1F("hChanSigCuts", "KKpi and piKK candidates",4,-0.5,3.5);
 362   fChanHistCuts[2] = new TH1F("hChanBkgCuts", "KKpi and piKK candidates",4,-0.5,3.5);
 363   fChanHistCuts[3] = new TH1F("hChanReflSigCuts", "KKpi and piKK candidates",4,-0.5,3.5);
 364   for(Int_t i=0;i<4;i++){
 365     fChanHist[i]->Sumw2();
 366     fChanHist[i]->SetMinimum(0);
 367     fChanHistCuts[i]->Sumw2();
 368     fChanHistCuts[i]->SetMinimum(0);
 369     fOutput->Add(fChanHist[i]);
 370     fOutput->Add(fChanHistCuts[i]);
 371   }
 372
 373   fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
 374   fHistNEvents->Sumw2();
 375   fHistNEvents->SetMinimum(0);
 376   fOutput->Add(fHistNEvents);
 377
 378
 379   return;
 380 }
 381
 382 //________________________________________________________________________
 383 void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
 384 {
 385   // Ds selection for current event, fill mass histos and selecetion variable histo
 386   // separate signal and backgound if fReadMC is activated
 387
 388   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 389   fHistNEvents->Fill(0); // count event
 390   // Post the data already here
 391   PostData(1,fOutput);
 392
 393
 394   TClonesArray *array3Prong = 0;
 395   if(!aod && AODEvent() && IsStandardAOD()) {
 396     // In case there is an AOD handler writing a standard AOD, use the AOD
 397     // event in memory rather than the input (ESD) event.
 398     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 399     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 400     // have to taken from the AOD event hold by the AliAODExtension
 401     AliAODHandler* aodHandler = (AliAODHandler*)
 402       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 403     if(aodHandler->GetExtensions()) {
 404       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 405       AliAODEvent *aodFromExt = ext->GetAOD();
 406       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 407     }
 408   } else {
 409     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 410   }
 411
 412   if(!array3Prong) {
 413     printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
 414     return;
 415   }
 416
 417
 418   TClonesArray *arrayMC=0;
 419   AliAODMCHeader *mcHeader=0;
 420
 421   // AOD primary vertex
 422   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 423   //    vtx1->Print();
 424
 425   // load MC particles
 426   if(fReadMC){
 427
 428     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 429     if(!arrayMC) {
 430       printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
 431       return;
 432     }
 433
 434   // load MC header
 435     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 436     if(!mcHeader) {
 437       printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
 438       return;
 439     }
 440   }
 441
 442   Int_t n3Prong = array3Prong->GetEntriesFast();
 443   if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
 444
 445
 446   Int_t pdgDstoKKpi[3]={321,321,211};
 447   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 448     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 449
 450
 451     Bool_t unsetvtx=kFALSE;
 452     if(!d->GetOwnPrimaryVtx()){
 453       d->SetOwnPrimaryVtx(vtx1);
 454       unsetvtx=kTRUE;
 455     }
 456
 457     Int_t retCodeProductionCuts=fProdCuts->IsSelected(d,AliRDHFCuts::kCandidate);
 458     if(retCodeProductionCuts>0){
 459       Int_t isKKpi=retCodeProductionCuts&1;
 460       Int_t ispiKK=retCodeProductionCuts&2;
 461 //       Int_t isPhi=retCodeProductionCuts&4;
 462 //       Int_t isK0star=retCodeProductionCuts&8;
 463       Double_t ptCand = d->Pt();
 464       Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
 465       Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate);
 466       Int_t isKKpiAC=retCodeAnalysisCuts&1;
 467       Int_t ispiKKAC=retCodeAnalysisCuts&2;
 468 //       Int_t isPhiAC=retCodeAnalysisCuts&4;
 469 //       Int_t isK0starAC=retCodeAnalysisCuts&8;
 470
 471       Int_t labDs=-1;
 472       if(fReadMC){
 473         labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
 474       }
 475
 476       Double_t dlen=d->DecayLength();
 477       Double_t cosp=d->CosPointingAngle();
 478       Double_t pt0=d->PtProng(0);
 479       Double_t pt1=d->PtProng(1);
 480       Double_t pt2=d->PtProng(2);
 481       Double_t sigvert=d->GetSigmaVert();
 482       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 483       Double_t dca=d->GetDCA();
 484       Double_t ptmax=0;
 485       for(Int_t i=0;i<3;i++){
 486         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 487       }
 488
 489       Int_t index=GetHistoIndex(iPtBin);
 490       Int_t type=0;
 491       if(isKKpi) type+=1;
 492       if(ispiKK) type+=2;
 493       Int_t typeAC=0;
 494       if(isKKpiAC) typeAC+=1;
 495       if(ispiKKAC) typeAC+=2;
 496       fCosPHist[index]->Fill(cosp);
 497       fDLenHist[index]->Fill(dlen);
 498       fSigVertHist[index]->Fill(sigvert);
 499       fSumd02Hist[index]->Fill(sumD02);
 500       fPtMaxHist[index]->Fill(ptmax);
 501       fPtCandHist[index]->Fill(ptCand);
 502       fDCAHist[index]->Fill(dca);
 503       fChanHist[0]->Fill(type);
 504       fPtProng0Hist[index]->Fill(pt0);
 505       fPtProng1Hist[index]->Fill(pt1);
 506       fPtProng2Hist[index]->Fill(pt2);
 507
 508       if(retCodeAnalysisCuts>0) fChanHistCuts[0]->Fill(typeAC);
 509       if(fReadMC){
 510         if(labDs>=0) {
 511           index=GetSignalHistoIndex(iPtBin);
 512           fChanHist[1]->Fill(type);
 513           if(retCodeAnalysisCuts>0) fChanHistCuts[1]->Fill(typeAC);
 514         }else{
 515           index=GetBackgroundHistoIndex(iPtBin);
 516           fChanHist[2]->Fill(type);
 517           if(retCodeAnalysisCuts>0) fChanHistCuts[2]->Fill(typeAC);
 518         }
 519       }
 520       if(isKKpi){
 521         index=GetHistoIndex(iPtBin);
 522         Double_t invMass=d->InvMassDsKKpi();
 523         fMassHist[index]->Fill(invMass);
 524         Double_t mass01=d->InvMass2Prongs(0,1,321,321);
 525         Double_t mass12=d->InvMass2Prongs(1,2,321,211);
 526         fDalitz[index]->Fill(mass01,mass12);
 527         if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
 528         if(fReadMC){
 529           Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 530           AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 531           Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 532           if(labDs>=0){
 533             if(pdgCode0==321) {
 534               index=GetSignalHistoIndex(iPtBin);
 535             }else{
 536               index=GetReflSignalHistoIndex(iPtBin);
 537             }
 538           }else{
 539             index=GetBackgroundHistoIndex(iPtBin);
 540           }
 541           fMassHist[index]->Fill(invMass);
 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           fPtCandHist[index]->Fill(ptCand);
 548           fDCAHist[index]->Fill(dca);
 549           fPtProng0Hist[index]->Fill(pt0);
 550           fPtProng1Hist[index]->Fill(pt1);
 551           fPtProng2Hist[index]->Fill(pt2);
 552           fDalitz[index]->Fill(mass01,mass12);
 553           if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
 554         }
 555       }
 556       if(ispiKK){
 557         index=GetHistoIndex(iPtBin);
 558         Double_t invMass=d->InvMassDspiKK();
 559         fMassHist[index]->Fill(invMass);
 560         if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
 561         if(fReadMC){
 562           Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 563           AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
 564           Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
 565           if(labDs>=0) {
 566             if(pdgCode0==211) {
 567               index=GetSignalHistoIndex(iPtBin);
 568             }else{
 569               index=GetReflSignalHistoIndex(iPtBin);
 570             }
 571           }else{
 572             index=GetBackgroundHistoIndex(iPtBin);
 573           }
 574           fMassHist[index]->Fill(invMass);
 575           fCosPHist[index]->Fill(cosp);
 576           fDLenHist[index]->Fill(dlen);
 577           fSigVertHist[index]->Fill(sigvert);
 578           fSumd02Hist[index]->Fill(sumD02);
 579           fPtMaxHist[index]->Fill(ptmax);
 580           fPtCandHist[index]->Fill(ptCand);
 581           fDCAHist[index]->Fill(dca);
 582           fPtProng0Hist[index]->Fill(pt0);
 583           fPtProng1Hist[index]->Fill(pt1);
 584           fPtProng2Hist[index]->Fill(pt2);
 585           if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
 586         }
 587       }
 588     }
 589     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 590   }
 591
 592
 593   PostData(1,fOutput);
 594   return;
 595 }
 596
 597 //_________________________________________________________________
 598
 599 void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
 600 {
 601   // Terminate analysis
 602   //
 603   if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
 604   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 605   if (!fOutput) {
 606     printf("ERROR: fOutput not available\n");
 607     return;
 608   }
 609   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
 610   fChanHist[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAll"));
 611   fChanHist[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSig"));
 612   fChanHist[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkg"));
 613   fChanHist[3] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanReflSig"));
 614   fChanHistCuts[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAllCuts"));
 615   fChanHistCuts[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSigCuts"));
 616   fChanHistCuts[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkgCuts"));
 617   fChanHistCuts[3] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanReflSigCuts"));
 618
 619
 620   TString hisname;
 621   Int_t index;
 622   for(Int_t i=0;i<fNPtBins;i++){
 623
 624     index=GetHistoIndex(i);
 625     hisname.Form("hMassAllPt%d",i);
 626     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 627     hisname.Form("hMassAllPt%dCuts",i);
 628     fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 629     hisname.Form("hCosPAllPt%d",i);
 630     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 631     hisname.Form("hDLenAllPt%d",i);
 632     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 633     hisname.Form("hSumd02AllPt%d",i);
 634     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 635     hisname.Form("hSigVertAllPt%d",i);
 636     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 637     hisname.Form("hPtMaxAllPt%d",i);
 638     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 639     hisname.Form("hPtCandAllPt%d",i);
 640     fPtCandHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 641     hisname.Form("hDCAAllPt%d",i);
 642     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 643     hisname.Form("hPtProng0AllPt%d",i);
 644     fPtProng0Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 645     hisname.Form("hPtProng1AllPt%d",i);
 646     fPtProng1Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 647     hisname.Form("hPtProng2AllPt%d",i);
 648     fPtProng2Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 649     hisname.Form("hDalitzAllPt%d",i);
 650     fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
 651
 652     index=GetSignalHistoIndex(i);
 653     hisname.Form("hMassSigPt%d",i);
 654     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 655     hisname.Form("hMassSigPt%dCuts",i);
 656     fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 657     hisname.Form("hCosPSigPt%d",i);
 658     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 659     hisname.Form("hDLenSigPt%d",i);
 660     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 661     hisname.Form("hSumd02SigPt%d",i);
 662     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 663     hisname.Form("hSigVertSigPt%d",i);
 664     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 665     hisname.Form("hPtMaxSigPt%d",i);
 666     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 667     hisname.Form("hPtCandSigPt%d",i);
 668     fPtCandHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 669     hisname.Form("hDCASigPt%d",i);
 670     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 671     hisname.Form("hPtProng0SigPt%d",i);
 672     fPtProng0Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 673     hisname.Form("hPtProng1SigPt%d",i);
 674     fPtProng1Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 675     hisname.Form("hPtProng2SigPt%d",i);
 676     fPtProng2Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 677     hisname.Form("hDalitzSigPt%d",i);
 678     fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
 679
 680     index=GetBackgroundHistoIndex(i);
 681     hisname.Form("hMassBkgPt%d",i);
 682     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 683     hisname.Form("hMassBkgPt%dCuts",i);
 684     fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 685     hisname.Form("hCosPBkgPt%d",i);
 686     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 687     hisname.Form("hDLenBkgPt%d",i);
 688     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 689     hisname.Form("hSumd02BkgPt%d",i);
 690     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 691     hisname.Form("hSigVertBkgPt%d",i);
 692     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 693     hisname.Form("hPtMaxBkgPt%d",i);
 694     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 695     hisname.Form("hPtCandBkgPt%d",i);
 696     fPtCandHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 697     hisname.Form("hDCABkgPt%d",i);
 698     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 699     hisname.Form("hPtProng0BkgPt%d",i);
 700     fPtProng0Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 701     hisname.Form("hPtProng1BkgPt%d",i);
 702     fPtProng1Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 703     hisname.Form("hPtProng2BkgPt%d",i);
 704     fPtProng2Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 705     hisname.Form("hDalitzBkgPt%d",i);
 706     fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
 707
 708     index=GetReflSignalHistoIndex(i);
 709     hisname.Form("hMassReflSigPt%d",i);
 710     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 711     hisname.Form("hMassReflSigPt%dCuts",i);
 712     fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 713     hisname.Form("hCosPReflSigPt%d",i);
 714     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 715     hisname.Form("hDLenReflSigPt%d",i);
 716     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 717     hisname.Form("hSumd02ReflSigPt%d",i);
 718     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 719     hisname.Form("hSigVertReflSigPt%d",i);
 720     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 721     hisname.Form("hPtMaxReflSigPt%d",i);
 722     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 723     hisname.Form("hPtCandReflSigPt%d",i);
 724     fPtCandHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 725     hisname.Form("hDCAReflSigPt%d",i);
 726     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 727     hisname.Form("hPtProng0ReflSigPt%d",i);
 728     fPtProng0Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 729     hisname.Form("hPtProng1ReflSigPt%d",i);
 730     fPtProng1Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 731     hisname.Form("hPtProng2ReflSigPt%d",i);
 732     fPtProng2Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 733     hisname.Form("hDalitzReflSigPt%d",i);
 734     fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
 735
 736   }
 737   return;
 738 }