PWGHF/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 "AliNormalizationCounter.h"
  45 #include "AliAnalysisTaskSEDs.h"
  46
  47 ClassImp(AliAnalysisTaskSEDs)
  48
  49
  50 //________________________________________________________________________
  51 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
  52   AliAnalysisTaskSE(),
  53   fOutput(0),
  54   fHistNEvents(0),
  55   fPtVsMass(0),
  56   fPtVsMassPhi(0),
  57   fPtVsMassK0st(0),
  58   fYVsPt(0),
  59   fYVsPtSig(0),
  60   fNtupleDs(0),
  61   fFillNtuple(0),
  62   fReadMC(kFALSE),
  63   fWriteOnlySignal(kFALSE),
  64   fDoCutVarHistos(kTRUE),
  65   fUseSelectionBit(kFALSE),
  66   fNPtBins(0),
  67   fListCuts(0),
  68   fMassRange(0.8),
  69   fMassBinSize(0.002),
  70   fCounter(0),
  71   fAnalysisCuts(0)
  72 {
  73   // Default constructor
  74
  75   for(Int_t i=0;i<3;i++){
  76     fHistCentrality[i]=0;
  77     fHistCentralityMult[i]=0;
  78   }
  79   for(Int_t i=0;i<4;i++) {
  80     fChanHist[i]=0;
  81   }
  82   for(Int_t i=0;i<4*kMaxPtBins;i++){
  83     fPtCandHist[i]=0;
  84     fMassHist[i]=0;
  85     fMassHistPhi[i]=0;
  86     fMassHistK0st[i]=0;
  87     fCosPHist[i]=0;
  88     fDLenHist[i]=0;
  89     fSumd02Hist[i]=0;
  90     fSigVertHist[i]=0;
  91     fPtMaxHist[i]=0;
  92     fDCAHist[i]=0;
  93     fPtProng0Hist[i]=0;
  94     fPtProng1Hist[i]=0;
  95     fPtProng2Hist[i]=0;
  96     fDalitz[i]=0;
  97     fDalitzPhi[i]=0;
  98     fDalitzK0st[i]=0;
  99   }
 100   for(Int_t i=0;i<kMaxPtBins;i++){
 101     fMassHistKK[i]=0;
 102     fMassHistKpi[i]=0;
 103   }
 104   for(Int_t i=0;i<kMaxPtBins+1;i++){
 105     fPtLimits[i]=0;
 106   }
 107
 108 }
 109
 110 //________________________________________________________________________
 111 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name,AliRDHFCutsDstoKKpi* analysiscuts,Int_t fillNtuple):
 112   AliAnalysisTaskSE(name),
 113   fOutput(0),
 114   fHistNEvents(0),
 115   fPtVsMass(0),
 116   fPtVsMassPhi(0),
 117   fPtVsMassK0st(0),
 118   fYVsPt(0),
 119   fYVsPtSig(0),
 120   fNtupleDs(0),
 121   fFillNtuple(fillNtuple),
 122   fReadMC(kFALSE),
 123   fWriteOnlySignal(kFALSE),
 124   fDoCutVarHistos(kTRUE),
 125   fUseSelectionBit(kFALSE),
 126   fNPtBins(0),
 127   fListCuts(0),
 128   fMassRange(0.8),
 129   fMassBinSize(0.002),
 130   fCounter(0),
 131   fAnalysisCuts(analysiscuts)
 132 {
 133   // Default constructor
 134   // Output slot #1 writes into a TList container
 135
 136   for(Int_t i=0;i<3;i++){
 137     fHistCentrality[i]=0;
 138     fHistCentralityMult[i]=0;
 139   }
 140   for(Int_t i=0;i<4;i++) {
 141     fChanHist[i]=0;
 142   }
 143   for(Int_t i=0;i<4*kMaxPtBins;i++){
 144     fPtCandHist[i]=0;
 145     fMassHist[i]=0;
 146     fMassHistPhi[i]=0;
 147     fMassHistK0st[i]=0;
 148     fCosPHist[i]=0;
 149     fDLenHist[i]=0;
 150     fSumd02Hist[i]=0;
 151     fSigVertHist[i]=0;
 152     fPtMaxHist[i]=0;
 153     fDCAHist[i]=0;
 154     fPtProng0Hist[i]=0;
 155     fPtProng1Hist[i]=0;
 156     fPtProng2Hist[i]=0;
 157     fDalitz[i]=0;
 158     fDalitzPhi[i]=0;
 159     fDalitzK0st[i]=0;
 160   }
 161   for(Int_t i=0;i<kMaxPtBins;i++){
 162     fMassHistKK[i]=0;
 163     fMassHistKpi[i]=0;
 164   }
 165   for(Int_t i=0;i<kMaxPtBins+1;i++){
 166     fPtLimits[i]=0;
 167   }
 168
 169   Int_t nptbins=fAnalysisCuts->GetNPtBins();
 170   Float_t *ptlim=fAnalysisCuts->GetPtBinLimits();
 171   SetPtBins(nptbins,ptlim);
 172
 173   DefineOutput(1,TList::Class());  //My private output
 174
 175   DefineOutput(2,TList::Class());
 176
 177   DefineOutput(3,AliNormalizationCounter::Class());
 178
 179   if(fFillNtuple>0){
 180     // Output slot #4 writes into a TNtuple container
 181     DefineOutput(4,TNtuple::Class());  //My private output
 182   }
 183
 184 }
 185
 186 //________________________________________________________________________
 187 void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Float_t* lim){
 188   // define pt bins for analysis
 189   if(n>kMaxPtBins){
 190     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
 191     fNPtBins=kMaxPtBins;
 192     fPtLimits[0]=0.;
 193     fPtLimits[1]=1.;
 194     fPtLimits[2]=3.;
 195     fPtLimits[3]=5.;
 196     fPtLimits[4]=10.;
 197     for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
 198   }else{
 199     fNPtBins=n;
 200     for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
 201     for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
 202   }
 203   if(fDebug > 1){
 204     printf("Number of Pt bins = %d\n",fNPtBins);
 205     for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);
 206   }
 207 }
 208 //________________________________________________________________________
 209 AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()
 210 {
 211   // Destructor
 212   if(fOutput && !fOutput->IsOwner()){
 213     delete fHistNEvents;
 214     for(Int_t i=0;i<4;i++){
 215       delete fChanHist[i];
 216     }
 217     for(Int_t i=0;i<4*fNPtBins;i++){
 218       delete fMassHist[i];
 219       delete fMassHistPhi[i];
 220       delete fMassHistK0st[i];
 221       delete fCosPHist[i];
 222       delete fDLenHist[i];
 223       delete fSumd02Hist[i];
 224       delete fSigVertHist[i];
 225       delete fPtMaxHist[i];
 226       delete fPtCandHist[i];
 227       delete fDCAHist[i];
 228       delete fPtProng0Hist[i];
 229       delete fPtProng1Hist[i];
 230       delete fPtProng2Hist[i];
 231       delete fDalitz[i];
 232       delete fDalitzPhi[i];
 233       delete fDalitzK0st[i];
 234     }
 235     for(Int_t i=0;i<fNPtBins;i++){
 236       delete fMassHistKK[i];
 237       delete fMassHistKpi[i];
 238     }
 239     delete fPtVsMass;
 240     delete fPtVsMassPhi;
 241     delete fPtVsMassK0st;
 242     delete fYVsPt;
 243     delete fYVsPtSig;
 244     for(Int_t i=0;i<3;i++){
 245       delete fHistCentrality[i];
 246       delete fHistCentralityMult[i];
 247     }
 248   }
 249   delete fOutput;
 250   delete fListCuts;
 251   delete fNtupleDs;
 252   delete fCounter;
 253   delete fAnalysisCuts;
 254
 255 }
 256
 257 //________________________________________________________________________
 258 void AliAnalysisTaskSEDs::Init()
 259 {
 260   // Initialization
 261
 262   if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
 263
 264   fListCuts=new TList();
 265   fListCuts->SetOwner();
 266   fListCuts->SetName("CutObjects");
 267
 268   AliRDHFCutsDstoKKpi *analysis = new AliRDHFCutsDstoKKpi(*fAnalysisCuts);
 269   analysis->SetName("AnalysisCuts");
 270
 271   fListCuts->Add(analysis);
 272   PostData(2,fListCuts);
 273   return;
 274 }
 275
 276 //________________________________________________________________________
 277 void AliAnalysisTaskSEDs::UserCreateOutputObjects()
 278 {
 279   // Create the output container
 280   //
 281   if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
 282
 283   // Several histograms are more conveniently managed in a TList
 284   fOutput = new TList();
 285   fOutput->SetOwner();
 286   fOutput->SetName("OutputHistos");
 287
 288   fHistNEvents = new TH1F("hNEvents", "number of events ",11,-0.5,10.5);
 289   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 290   fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
 291   fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
 292   fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to not reco vertex");
 293   fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected for contr vertex");
 294   fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept");
 295   fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events");
 296   fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
 297   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of 3 prong candidates");
 298   fHistNEvents->GetXaxis()->SetBinLabel(10,"no. of Ds after filtering cuts");
 299   fHistNEvents->GetXaxis()->SetBinLabel(11,"no. of Ds after selection cuts");
 300
 301   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 302
 303   fHistNEvents->Sumw2();
 304   fHistNEvents->SetMinimum(0);
 305   fOutput->Add(fHistNEvents);
 306
 307   fHistCentrality[0]=new TH1F("hCentr","centrality",10000,0.,100.);
 308   fHistCentrality[1]=new TH1F("hCentr(selectedCent)","centrality(selectedCent)",10000,0.,100.);
 309   fHistCentrality[2]=new TH1F("hCentr(OutofCent)","centrality(OutofCent)",10000,0.,100.);
 310   fHistCentralityMult[0]=new TH2F("hCentrMult","centrality vs mult",100,0.5,30000.5,40,0.,100.);
 311   fHistCentralityMult[1]=new TH2F("hCentrMult(selectedCent)","centrality vs mult(selectedCent)",100,0.5,30000.5,40,0.,100.);
 312   fHistCentralityMult[2]=new TH2F("hCentrMult(OutofCent)","centrality vs mult(OutofCent)",100,0.5,30000.5,40,0.,100.);
 313   for(Int_t i=0;i<3;i++){
 314     fHistCentrality[i]->Sumw2();
 315     fOutput->Add(fHistCentrality[i]);
 316     fHistCentralityMult[i]->Sumw2();
 317     fOutput->Add(fHistCentralityMult[i]);
 318   }
 319
 320   Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
 321
 322   Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
 323   if(nInvMassBins%2==1) nInvMassBins++;
 324   // Double_t minMass=massDs-1.0*nInvMassBins*fMassBinSize;
 325   Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
 326   //  Double_t maxMass=massDs+1.0*nInvMassBins*fMassBinSize;
 327   Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
 328
 329   TString hisname;
 330   TString htype;
 331   Int_t index;
 332   for(Int_t iType=0; iType<4; iType++){
 333     for(Int_t i=0;i<fNPtBins;i++){
 334       if(iType==0){
 335         htype="All";
 336         index=GetHistoIndex(i);
 337       }else if(iType==1){
 338         htype="Sig";
 339         index=GetSignalHistoIndex(i);
 340       }else if(iType==2){
 341         htype="Bkg";
 342         index=GetBackgroundHistoIndex(i);
 343       }else{
 344         htype="ReflSig";
 345         index=GetReflSignalHistoIndex(i);
 346       }
 347       hisname.Form("hMass%sPt%d",htype.Data(),i);
 348       fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 349       fMassHist[index]->Sumw2();
 350       hisname.Form("hMass%sPt%dphi",htype.Data(),i);
 351       fMassHistPhi[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 352       fMassHistPhi[index]->Sumw2();
 353       hisname.Form("hMass%sPt%dk0st",htype.Data(),i);
 354       fMassHistK0st[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
 355       fMassHistK0st[index]->Sumw2();
 356       hisname.Form("hCosP%sPt%d",htype.Data(),i);
 357       fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 358       fCosPHist[index]->Sumw2();
 359       hisname.Form("hDLen%sPt%d",htype.Data(),i);
 360       fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 361       fDLenHist[index]->Sumw2();
 362       hisname.Form("hSumd02%sPt%d",htype.Data(),i);
 363       fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 364       fSumd02Hist[index]->Sumw2();
 365       hisname.Form("hSigVert%sPt%d",htype.Data(),i);
 366       fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 367       fSigVertHist[index]->Sumw2();
 368       hisname.Form("hPtMax%sPt%d",htype.Data(),i);
 369       fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
 370       fPtMaxHist[index]->Sumw2();
 371       hisname.Form("hPtCand%sPt%d",htype.Data(),i);
 372       fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
 373       fPtCandHist[index]->Sumw2();
 374       hisname.Form("hDCA%sPt%d",htype.Data(),i);
 375       fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 376       fDCAHist[index]->Sumw2();
 377       hisname.Form("hPtProng0%sPt%d",htype.Data(),i);
 378       fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 379       fPtProng0Hist[index]->Sumw2();
 380       hisname.Form("hPtProng1%sPt%d",htype.Data(),i);
 381       fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 382       fPtProng1Hist[index]->Sumw2();
 383       hisname.Form("hPtProng2%sPt%d",htype.Data(),i);
 384       fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
 385       fPtProng2Hist[index]->Sumw2();
 386       hisname.Form("hDalitz%sPt%d",htype.Data(),i);
 387       fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 388       fDalitz[index]->Sumw2();
 389       hisname.Form("hDalitz%sPt%dphi",htype.Data(),i);
 390       fDalitzPhi[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 391       fDalitzPhi[index]->Sumw2();
 392       hisname.Form("hDalitz%sPt%dk0st",htype.Data(),i);
 393       fDalitzK0st[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 394       fDalitzK0st[index]->Sumw2();
 395     }
 396   }
 397
 398   for(Int_t i=0; i<4*fNPtBins; i++){
 399     fOutput->Add(fMassHist[i]);
 400     fOutput->Add(fMassHistPhi[i]);
 401     fOutput->Add(fMassHistK0st[i]);
 402     fOutput->Add(fPtCandHist[i]);
 403     if(fDoCutVarHistos){
 404       fOutput->Add(fCosPHist[i]);
 405       fOutput->Add(fDLenHist[i]);
 406       fOutput->Add(fSumd02Hist[i]);
 407       fOutput->Add(fSigVertHist[i]);
 408       fOutput->Add(fPtMaxHist[i]);
 409       fOutput->Add(fDCAHist[i]);
 410       fOutput->Add(fPtProng0Hist[i]);
 411       fOutput->Add(fPtProng1Hist[i]);
 412       fOutput->Add(fPtProng2Hist[i]);
 413       fOutput->Add(fDalitz[i]);
 414       fOutput->Add(fDalitzPhi[i]);
 415       fOutput->Add(fDalitzK0st[i]);
 416     }
 417   }
 418
 419   fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",64,-0.5,63.5);
 420   fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",64,-0.5,63.5);
 421   fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",64,-0.5,63.5);
 422   fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",64,-0.5,63.5);
 423   for(Int_t i=0;i<4;i++){
 424     fChanHist[i]->Sumw2();
 425     fChanHist[i]->SetMinimum(0);
 426     fOutput->Add(fChanHist[i]);
 427   }
 428
 429
 430   fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nInvMassBins,minMass,maxMass,40,0.,20.);
 431   fPtVsMassPhi=new TH2F("hPtVsMassPhi","PtVsMass (phi selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
 432   fPtVsMassK0st=new TH2F("hPtVsMassK0st","PtVsMass (K0* selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
 433   fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
 434   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
 435
 436   for(Int_t i=0;i<fNPtBins;i++){
 437     hisname.Form("hMassKKPt%d",i);
 438     fMassHistKK[i]=new TH1F(hisname.Data(),hisname.Data(),200,0.95,1.35);
 439     fMassHistKK[i]->Sumw2();
 440     fOutput->Add(fMassHistKK[i]);
 441     hisname.Form("hMassKpiPt%d",i);
 442     fMassHistKpi[i]=new TH1F(hisname.Data(),hisname.Data(),200,0.7,1.1);
 443     fMassHistKpi[i]->Sumw2();
 444     fOutput->Add(fMassHistKpi[i]);
 445   }
 446
 447   fOutput->Add(fPtVsMass);
 448   fOutput->Add(fPtVsMassPhi);
 449   fOutput->Add(fPtVsMassK0st);
 450   fOutput->Add(fYVsPt);
 451   fOutput->Add(fYVsPtSig);
 452
 453   //Counter for Normalization
 454   fCounter = new AliNormalizationCounter("NormalizationCounter");
 455   fCounter->Init();
 456
 457   PostData(1,fOutput);
 458   PostData(3,fCounter);
 459
 460   if(fFillNtuple>0){
 461     OpenFile(4); // 4 is the slot number of the ntuple
 462
 463     fNtupleDs = new TNtuple("fNtupleDs","Ds","labDs:retcode:pdgcode0:Pt0:Pt1:Pt2:PtRec:P0:P1:P2:PidTrackBit0:PidTrackBit1:PidTrackBit2:PointingAngle:PointingAngleXY:DecLeng:DecLengXY:NorDecLeng:NorDecLengXY:InvMassKKpi:InvMasspiKK:sigvert:d00:d01:d02:dca:d0square:InvMassPhiKKpi:InvMassPhipiKK:InvMassK0starKKpi:InvMassK0starpiKK:cosinePiDsFrameKKpi:cosinePiDsFramepiKK:cosineKPhiFrameKKpi:cosineKPhiFramepiKK:centrality:runNumber");
 464
 465   }
 466
 467
 468   return;
 469 }
 470
 471 //________________________________________________________________________
 472 void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
 473 {
 474   // Ds selection for current event, fill mass histos and selecetion variable histo
 475   // separate signal and backgound if fReadMC is activated
 476
 477   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 478
 479   TClonesArray *array3Prong = 0;
 480   if(!aod && AODEvent() && IsStandardAOD()) {
 481     // In case there is an AOD handler writing a standard AOD, use the AOD
 482     // event in memory rather than the input (ESD) event.
 483     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 484     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 485     // have to taken from the AOD event hold by the AliAODExtension
 486     AliAODHandler* aodHandler = (AliAODHandler*)
 487       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 488     if(aodHandler->GetExtensions()) {
 489       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 490       AliAODEvent *aodFromExt = ext->GetAOD();
 491       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 492     }
 493   } else if(aod) {
 494     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 495   }
 496
 497   if(!aod || !array3Prong) {
 498     printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
 499     return;
 500   }
 501
 502
 503   // fix for temporary bug in ESDfilter
 504   // the AODs with null vertex pointer didn't pass the PhysSel
 505   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 506
 507
 508   fHistNEvents->Fill(0); // count event
 509   // Post the data already here
 510   PostData(1,fOutput);
 511
 512   fCounter->StoreEvent(aod,fAnalysisCuts,fReadMC);
 513
 514
 515   Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
 516   Float_t ntracks=aod->GetNTracks();
 517   Float_t evCentr=fAnalysisCuts->GetCentrality(aod);
 518
 519   fHistCentrality[0]->Fill(evCentr);
 520   fHistCentralityMult[0]->Fill(ntracks,evCentr);
 521   if(fAnalysisCuts->IsEventRejectedDueToTrigger())fHistNEvents->Fill(2);
 522   if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3);
 523   if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4);
 524   if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5);
 525   if(fAnalysisCuts->IsEventRejectedDueToPileup())fHistNEvents->Fill(6);
 526   if(fAnalysisCuts->IsEventRejectedDueToCentrality()){
 527     fHistNEvents->Fill(7);
 528     fHistCentrality[2]->Fill(evCentr);
 529     fHistCentralityMult[2]->Fill(ntracks,evCentr);
 530   }
 531
 532   Float_t centrality=fAnalysisCuts->GetCentrality(aod);
 533   Int_t runNumber=aod->GetRunNumber();
 534
 535   if(!isEvSel)return;
 536
 537   fHistNEvents->Fill(1);
 538   fHistCentrality[1]->Fill(evCentr);
 539   fHistCentralityMult[1]->Fill(ntracks,evCentr);
 540
 541   TClonesArray *arrayMC=0;
 542   AliAODMCHeader *mcHeader=0;
 543
 544   // AOD primary vertex
 545   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 546   //    vtx1->Print();
 547
 548   // load MC particles
 549   if(fReadMC){
 550
 551     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 552     if(!arrayMC) {
 553       printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
 554       return;
 555     }
 556
 557   // load MC header
 558     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 559     if(!mcHeader) {
 560       printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
 561       return;
 562     }
 563   }
 564
 565   Int_t n3Prong = array3Prong->GetEntriesFast();
 566   if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
 567
 568
 569   Int_t pdgDstoKKpi[3]={321,321,211};
 570   Int_t nSelected=0;
 571   Int_t nFiltered=0;
 572   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 573
 574     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 575     fHistNEvents->Fill(8);
 576
 577     if(fUseSelectionBit && !(d->HasSelectionBit(AliRDHFCuts::kDsCuts))){
 578       continue;
 579     }
 580     nFiltered++;
 581     fHistNEvents->Fill(9);
 582
 583     Bool_t unsetvtx=kFALSE;
 584     if(!d->GetOwnPrimaryVtx()){
 585       d->SetOwnPrimaryVtx(vtx1);
 586       unsetvtx=kTRUE;
 587     }
 588
 589     Bool_t recVtx=kFALSE;
 590     AliAODVertex *origownvtx=0x0;
 591
 592
 593     Double_t ptCand = d->Pt();
 594     Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
 595     Double_t rapid=d->YDs();
 596     fYVsPt->Fill(ptCand,rapid);
 597     Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
 598     if(!isFidAcc) continue;
 599
 600     Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
 601    Int_t retCodeNoRes=retCodeAnalysisCuts;
 602     Bool_t origRes=fAnalysisCuts->IsCutOnResonancesApplied();
 603     if(origRes){
 604       fAnalysisCuts->ApplyCutOnResonances(kFALSE);
 605       retCodeNoRes=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
 606       fAnalysisCuts->ApplyCutOnResonances(origRes);
 607     }
 608     if(retCodeNoRes&1){ //KKpi
 609       Double_t massKK=d->InvMass2Prongs(0,1,321,321);
 610       Double_t massKp=d->InvMass2Prongs(1,2,321,211);
 611       fMassHistKK[iPtBin]->Fill(massKK);
 612       fMassHistKpi[iPtBin]->Fill(massKp);
 613     }
 614     if(retCodeNoRes&2){ //piKK
 615       Double_t massKK=d->InvMass2Prongs(1,2,321,321);
 616       Double_t massKp=d->InvMass2Prongs(0,1,211,321);
 617       fMassHistKK[iPtBin]->Fill(massKK);
 618       fMassHistKpi[iPtBin]->Fill(massKp);
 619     }
 620
 621     if(retCodeAnalysisCuts<=0) continue;
 622
 623     if(fAnalysisCuts->GetIsPrimaryWithoutDaughters()){
 624       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
 625       if(fAnalysisCuts->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
 626       else fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
 627     }
 628
 629
 630     fHistNEvents->Fill(10);
 631     nSelected++;
 632
 633     Int_t index=GetHistoIndex(iPtBin);
 634     fPtCandHist[index]->Fill(ptCand);
 635
 636     Int_t isKKpi=retCodeAnalysisCuts&1;
 637     Int_t ispiKK=retCodeAnalysisCuts&2;
 638     Int_t isPhiKKpi=retCodeAnalysisCuts&4;
 639     Int_t isPhipiKK=retCodeAnalysisCuts&8;
 640     Int_t isK0starKKpi=retCodeAnalysisCuts&16;
 641     Int_t isK0starpiKK=retCodeAnalysisCuts&32;
 642
 643     Double_t weightKKpi=1.;
 644     Double_t weightpiKK=1.;
 645     if(fAnalysisCuts->GetPidOption()==AliRDHFCutsDstoKKpi::kBayesianWeights){
 646       weightKKpi=fAnalysisCuts->GetWeightForKKpi();
 647       weightpiKK=fAnalysisCuts->GetWeightForpiKK();
 648       if(weightKKpi>1. || weightKKpi<0.) weightKKpi=0.;
 649       if(weightpiKK>1. || weightpiKK<0.) weightpiKK=0.;
 650     }
 651
 652     fChanHist[0]->Fill(retCodeAnalysisCuts);
 653
 654     Int_t indexMCKKpi=-1;
 655     Int_t indexMCpiKK=-1;
 656     Int_t labDs=-1;
 657     Int_t pdgCode0=-999;
 658     Int_t isMCSignal=-1;
 659
 660     if(fReadMC){
 661       labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
 662       if(labDs>=0){
 663         Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 664         AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
 665         pdgCode0=TMath::Abs(p->GetPdgCode());
 666
 667         if(isKKpi){
 668           if(pdgCode0==321) {
 669             indexMCKKpi=GetSignalHistoIndex(iPtBin);
 670             fYVsPtSig->Fill(ptCand,rapid);
 671             fChanHist[1]->Fill(retCodeAnalysisCuts);
 672             isMCSignal=1;
 673           }else{
 674             indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
 675             fChanHist[3]->Fill(retCodeAnalysisCuts);
 676             isMCSignal=0;
 677           }
 678         }
 679         if(ispiKK){
 680           if(pdgCode0==211) {
 681             indexMCpiKK=GetSignalHistoIndex(iPtBin);
 682             fYVsPtSig->Fill(ptCand,rapid);
 683             fChanHist[1]->Fill(retCodeAnalysisCuts);
 684             isMCSignal=1;
 685           }else{
 686             indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
 687             fChanHist[3]->Fill(retCodeAnalysisCuts);
 688             isMCSignal=0;
 689           }
 690         }
 691       }else{
 692         indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
 693         indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
 694         fChanHist[2]->Fill(retCodeAnalysisCuts);
 695       }
 696     }
 697
 698     if(isKKpi){
 699       Double_t invMass=d->InvMassDsKKpi();
 700       fMassHist[index]->Fill(invMass,weightKKpi);
 701       fPtVsMass->Fill(invMass,ptCand,weightKKpi);
 702       if(isPhiKKpi){
 703         fMassHistPhi[index]->Fill(invMass,weightKKpi);
 704         fPtVsMassPhi->Fill(invMass,ptCand,weightKKpi);
 705       }
 706       if(isK0starKKpi){
 707         fMassHistK0st[index]->Fill(invMass,weightKKpi);
 708         fPtVsMassK0st->Fill(invMass,ptCand,weightKKpi);
 709       }
 710       if(fReadMC  && indexMCKKpi!=-1){
 711         fMassHist[indexMCKKpi]->Fill(invMass,weightKKpi);
 712         if(isPhiKKpi) fMassHistPhi[indexMCKKpi]->Fill(invMass,weightKKpi);
 713         if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass,weightKKpi);
 714       }
 715     }
 716     if(ispiKK){
 717       Double_t invMass=d->InvMassDspiKK();
 718       fMassHist[index]->Fill(invMass,weightpiKK);
 719       fPtVsMass->Fill(invMass,ptCand,weightpiKK);
 720       if(isPhipiKK){
 721         fMassHistPhi[index]->Fill(invMass,weightpiKK);
 722         fPtVsMassPhi->Fill(invMass,ptCand,weightpiKK);
 723       }
 724       if(isK0starpiKK){
 725         fMassHistK0st[index]->Fill(invMass,weightpiKK);
 726         fPtVsMassK0st->Fill(invMass,ptCand,weightpiKK);
 727       }
 728       if(fReadMC  && indexMCpiKK!=-1){
 729         fMassHist[indexMCpiKK]->Fill(invMass,weightpiKK);
 730         if(isPhipiKK) fMassHistPhi[indexMCpiKK]->Fill(invMass,weightpiKK);
 731         if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass,weightpiKK);
 732       }
 733     }
 734
 735     if(fDoCutVarHistos){
 736       Double_t dlen=d->DecayLength();
 737       Double_t cosp=d->CosPointingAngle();
 738       Double_t pt0=d->PtProng(0);
 739       Double_t pt1=d->PtProng(1);
 740       Double_t pt2=d->PtProng(2);
 741       Double_t sigvert=d->GetSigmaVert();
 742       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 743       Double_t dca=d->GetDCA();
 744       Double_t ptmax=0;
 745       for(Int_t i=0;i<3;i++){
 746         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 747       }
 748       fCosPHist[index]->Fill(cosp);
 749       fDLenHist[index]->Fill(dlen);
 750       fSigVertHist[index]->Fill(sigvert);
 751       fSumd02Hist[index]->Fill(sumD02);
 752       fPtMaxHist[index]->Fill(ptmax);
 753       fDCAHist[index]->Fill(dca);
 754       fPtProng0Hist[index]->Fill(pt0);
 755       fPtProng1Hist[index]->Fill(pt1);
 756       fPtProng2Hist[index]->Fill(pt2);
 757       if(isKKpi){
 758         Double_t massKK=d->InvMass2Prongs(0,1,321,321);
 759         Double_t massKp=d->InvMass2Prongs(1,2,321,211);
 760         fDalitz[index]->Fill(massKK,massKp);
 761         if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
 762         if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
 763         if(fReadMC && indexMCKKpi!=-1){
 764           fDalitz[indexMCKKpi]->Fill(massKK,massKp);
 765           if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
 766           if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
 767           fCosPHist[indexMCKKpi]->Fill(cosp);
 768           fDLenHist[indexMCKKpi]->Fill(dlen);
 769           fSigVertHist[indexMCKKpi]->Fill(sigvert);
 770           fSumd02Hist[indexMCKKpi]->Fill(sumD02);
 771           fPtMaxHist[indexMCKKpi]->Fill(ptmax);
 772           fPtCandHist[indexMCKKpi]->Fill(ptCand);
 773           fDCAHist[indexMCKKpi]->Fill(dca);
 774           fPtProng0Hist[indexMCKKpi]->Fill(pt0);
 775           fPtProng1Hist[indexMCKKpi]->Fill(pt1);
 776           fPtProng2Hist[indexMCKKpi]->Fill(pt2);
 777         }
 778       }
 779       if(ispiKK){
 780         Double_t massKK=d->InvMass2Prongs(1,2,321,321);
 781         Double_t massKp=d->InvMass2Prongs(0,1,211,321);
 782         fDalitz[index]->Fill(massKK,massKp);
 783         if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
 784         if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
 785
 786
 787         if(fReadMC && indexMCpiKK!=-1){
 788           fDalitz[indexMCpiKK]->Fill(massKK,massKp);
 789           if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
 790           if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
 791           fCosPHist[indexMCpiKK]->Fill(cosp);
 792           fDLenHist[indexMCpiKK]->Fill(dlen);
 793           fSigVertHist[indexMCpiKK]->Fill(sigvert);
 794           fSumd02Hist[indexMCpiKK]->Fill(sumD02);
 795           fPtMaxHist[indexMCpiKK]->Fill(ptmax);
 796           fPtCandHist[indexMCpiKK]->Fill(ptCand);
 797           fDCAHist[indexMCpiKK]->Fill(dca);
 798           fPtProng0Hist[indexMCpiKK]->Fill(pt0);
 799           fPtProng1Hist[indexMCpiKK]->Fill(pt1);
 800           fPtProng2Hist[indexMCpiKK]->Fill(pt2);
 801         }
 802       }
 803
 804
 805     }
 806
 807     Float_t tmp[37];
 808     if(fFillNtuple>0){
 809
 810       if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){
 811
 812         AliAODTrack *track0=(AliAODTrack*)d->GetDaughter(0);
 813         AliAODTrack *track1=(AliAODTrack*)d->GetDaughter(1);
 814         AliAODTrack *track2=(AliAODTrack*)d->GetDaughter(2);
 815
 816         UInt_t BitMapPIDTrack0=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track0);
 817         UInt_t BitMapPIDTrack1=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track1);
 818         UInt_t BitMapPIDTrack2=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track2);
 819
 820             tmp[0]=Float_t(labDs);
 821             if(fReadMC && fWriteOnlySignal) tmp[0]=Float_t(isMCSignal);
 822             tmp[1]=Float_t(retCodeAnalysisCuts);
 823             tmp[2]=Float_t(pdgCode0);
 824             tmp[3]=d->PtProng(0);
 825             tmp[4]=d->PtProng(1);
 826             tmp[5]=d->PtProng(2);
 827             tmp[6]=d->Pt();
 828             tmp[7]=d->PProng(0);
 829             tmp[8]=d->PProng(1);
 830             tmp[9]=d->PProng(2);
 831             tmp[10]=Int_t(BitMapPIDTrack0);
 832             tmp[11]=Int_t(BitMapPIDTrack1);
 833             tmp[12]=Int_t(BitMapPIDTrack2);
 834             tmp[13]=d->CosPointingAngle();
 835             tmp[14]=d->CosPointingAngleXY();
 836             tmp[15]=d->DecayLength();
 837             tmp[16]=d->DecayLengthXY();
 838             tmp[17]=d->NormalizedDecayLength();
 839             tmp[18]=d->NormalizedDecayLengthXY();
 840             tmp[19]=d->InvMassDsKKpi();
 841             tmp[20]=d->InvMassDspiKK();
 842             tmp[21]=d->GetSigmaVert();
 843             tmp[22]=d->Getd0Prong(0);
 844             tmp[23]=d->Getd0Prong(1);
 845             tmp[24]=d->Getd0Prong(2);
 846             tmp[25]=d->GetDCA();
 847             tmp[26]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 848             tmp[27]=d->InvMass2Prongs(0,1,321,321);
 849             tmp[28]=d->InvMass2Prongs(1,2,321,321);
 850             tmp[29]=d->InvMass2Prongs(1,2,321,211);
 851             tmp[30]=d->InvMass2Prongs(0,1,211,321);
 852             tmp[31]=d->CosPiDsLabFrameKKpi();
 853             tmp[32]=d->CosPiDsLabFramepiKK();
 854             tmp[33]=d->CosPiKPhiRFrameKKpi();
 855             tmp[34]=d->CosPiKPhiRFramepiKK();
 856             tmp[35]=(Float_t)(centrality);
 857             tmp[36]=(Float_t)(runNumber);
 858
 859             if(fReadMC && fWriteOnlySignal){
 860               if(isMCSignal>=0) fNtupleDs->Fill(tmp);
 861             }else{
 862               fNtupleDs->Fill(tmp);
 863             }
 864             PostData(4,fNtupleDs);
 865       }
 866     }
 867
 868     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 869     if(recVtx)fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
 870   }
 871
 872   fCounter->StoreCandidates(aod,nFiltered,kTRUE);
 873   fCounter->StoreCandidates(aod,nSelected,kFALSE);
 874
 875   PostData(1,fOutput);
 876   PostData(3,fCounter);
 877
 878   return;
 879 }
 880
 881 //_________________________________________________________________
 882
 883 void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
 884 {
 885   // Terminate analysis
 886   //
 887   if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
 888   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 889   if (!fOutput) {
 890     printf("ERROR: fOutput not available\n");
 891     return;
 892   }
 893   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
 894   if(fHistNEvents){
 895     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
 896   }else{
 897     printf("ERROR: fHistNEvents not available\n");
 898     return;
 899   }
 900   return;
 901 }
 902