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