PWG3/vertexingHF/AliAnalysisTaskSEDs.cxx

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