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.4),
  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.4),
  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
 378   PostData(1,fOutput);
 379   PostData(3,fCounter);
 380
 381   if(fFillNtuple>0){
 382     OpenFile(4); // 4 is the slot number of the ntuple
 383
 384     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");
 385
 386   }
 387
 388
 389   return;
 390 }
 391
 392 //________________________________________________________________________
 393 void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
 394 {
 395   // Ds selection for current event, fill mass histos and selecetion variable histo
 396   // separate signal and backgound if fReadMC is activated
 397
 398   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 399
 400   TClonesArray *array3Prong = 0;
 401   if(!aod && AODEvent() && IsStandardAOD()) {
 402     // In case there is an AOD handler writing a standard AOD, use the AOD
 403     // event in memory rather than the input (ESD) event.
 404     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 405     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 406     // have to taken from the AOD event hold by the AliAODExtension
 407     AliAODHandler* aodHandler = (AliAODHandler*)
 408       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 409     if(aodHandler->GetExtensions()) {
 410       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 411       AliAODEvent *aodFromExt = ext->GetAOD();
 412       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 413     }
 414   } else if(aod) {
 415     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 416   }
 417
 418   if(!aod || !array3Prong) {
 419     printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
 420     return;
 421   }
 422
 423
 424   // fix for temporary bug in ESDfilter
 425   // the AODs with null vertex pointer didn't pass the PhysSel
 426   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 427
 428
 429   fHistNEvents->Fill(0); // count event
 430   // Post the data already here
 431   PostData(1,fOutput);
 432
 433   fCounter->StoreEvent(aod,fProdCuts,fReadMC);
 434
 435
 436   Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
 437   if(fAnalysisCuts->IsEventRejectedDueToTrigger())fHistNEvents->Fill(2);
 438   if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3);
 439   if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4);
 440   if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5);
 441   if(fAnalysisCuts->IsEventRejectedDueToPileupSPD())fHistNEvents->Fill(6);
 442   if(fAnalysisCuts->IsEventRejectedDueToCentrality())fHistNEvents->Fill(7);
 443
 444
 445
 446   if(!isEvSel)return;
 447
 448   fHistNEvents->Fill(1);
 449
 450   TClonesArray *arrayMC=0;
 451   AliAODMCHeader *mcHeader=0;
 452
 453   // AOD primary vertex
 454   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 455   //    vtx1->Print();
 456
 457   // load MC particles
 458   if(fReadMC){
 459
 460     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 461     if(!arrayMC) {
 462       printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
 463       return;
 464     }
 465
 466   // load MC header
 467     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 468     if(!mcHeader) {
 469       printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
 470       return;
 471     }
 472   }
 473
 474   Int_t n3Prong = array3Prong->GetEntriesFast();
 475   if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
 476
 477
 478   Int_t pdgDstoKKpi[3]={321,321,211};
 479   Int_t nSelectedloose=0;
 480   Int_t nSelectedtight=0;
 481
 482   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 483
 484     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 485     fHistNEvents->Fill(8);
 486
 487     if(fUseSelectionBit && !(d->HasSelectionBit(AliRDHFCuts::kDsCuts))){
 488       fHistNEvents->Fill(11);
 489       continue;
 490     }
 491
 492     Bool_t unsetvtx=kFALSE;
 493     if(!d->GetOwnPrimaryVtx()){
 494       d->SetOwnPrimaryVtx(vtx1);
 495       unsetvtx=kTRUE;
 496     }
 497
 498     Bool_t recVtx=kFALSE;
 499     AliAODVertex *origownvtx=0x0;
 500     Int_t retCodeProdCuts=fProdCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 501
 502     if(retCodeProdCuts) {
 503       if(fProdCuts->GetIsPrimaryWithoutDaughters()){
 504             if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
 505             if(fProdCuts->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
 506             else fProdCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
 507       }
 508     }
 509
 510     Double_t ptCand = d->Pt();
 511     Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
 512     Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 513     Double_t rapid=d->YDs();
 514     fYVsPt->Fill(ptCand,rapid);
 515
 516     Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
 517
 518     if(retCodeProdCuts>0){
 519       if(isFidAcc){
 520         nSelectedloose++;
 521         fHistNEvents->Fill(9);
 522         if(retCodeAnalysisCuts>0)nSelectedtight++;
 523       }
 524     }
 525
 526     if(retCodeAnalysisCuts<=0) continue;
 527     if(!isFidAcc) continue;
 528     fHistNEvents->Fill(10);
 529
 530     Int_t index=GetHistoIndex(iPtBin);
 531     fPtCandHist[index]->Fill(ptCand);
 532
 533     Int_t isKKpi=retCodeAnalysisCuts&1;
 534     Int_t ispiKK=retCodeAnalysisCuts&2;
 535     Int_t isPhiKKpi=retCodeAnalysisCuts&4;
 536     Int_t isPhipiKK=retCodeAnalysisCuts&8;
 537     Int_t isK0starKKpi=retCodeAnalysisCuts&16;
 538     Int_t isK0starpiKK=retCodeAnalysisCuts&32;
 539
 540     fChanHist[0]->Fill(retCodeAnalysisCuts);
 541
 542     Int_t indexMCKKpi=-1;
 543     Int_t indexMCpiKK=-1;
 544     Int_t labDs=-1;
 545     Int_t pdgCode0=-999;
 546
 547     if(fReadMC){
 548       labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
 549       if(labDs>=0){
 550         Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
 551         AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
 552         pdgCode0=TMath::Abs(p->GetPdgCode());
 553
 554         if(isKKpi){
 555           if(pdgCode0==321) {
 556             indexMCKKpi=GetSignalHistoIndex(iPtBin);
 557             fYVsPtSig->Fill(ptCand,rapid);
 558             fChanHist[1]->Fill(retCodeAnalysisCuts);
 559           }else{
 560             indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
 561             fChanHist[3]->Fill(retCodeAnalysisCuts);
 562           }
 563         }
 564         if(ispiKK){
 565           if(pdgCode0==211) {
 566             indexMCpiKK=GetSignalHistoIndex(iPtBin);
 567             fYVsPtSig->Fill(ptCand,rapid);
 568             fChanHist[1]->Fill(retCodeAnalysisCuts);
 569           }else{
 570             indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
 571             fChanHist[3]->Fill(retCodeAnalysisCuts);
 572           }
 573         }
 574       }else{
 575         indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
 576         indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
 577         fChanHist[2]->Fill(retCodeAnalysisCuts);
 578       }
 579     }
 580
 581     if(isKKpi){
 582       Double_t invMass=d->InvMassDsKKpi();
 583       fMassHist[index]->Fill(invMass);
 584       fPtVsMass->Fill(invMass,ptCand);
 585       if(isPhiKKpi) fMassHistPhi[index]->Fill(invMass);
 586       if(isK0starKKpi) fMassHistK0st[index]->Fill(invMass);
 587       if(fReadMC  && indexMCKKpi!=-1){
 588         fMassHist[indexMCKKpi]->Fill(invMass);
 589         if(isPhiKKpi) fMassHistPhi[indexMCKKpi]->Fill(invMass);
 590         if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass);
 591       }
 592     }
 593     if(ispiKK){
 594       Double_t invMass=d->InvMassDspiKK();
 595       fMassHist[index]->Fill(invMass);
 596       fPtVsMass->Fill(invMass,ptCand);
 597       if(isPhipiKK) fMassHistPhi[index]->Fill(invMass);
 598       if(isK0starpiKK) fMassHistK0st[index]->Fill(invMass);
 599       if(fReadMC  && indexMCpiKK!=-1){
 600         fMassHist[indexMCpiKK]->Fill(invMass);
 601         if(isPhipiKK) fMassHistPhi[indexMCpiKK]->Fill(invMass);
 602         if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass);
 603       }
 604     }
 605
 606     if(fDoCutVarHistos){
 607       Double_t dlen=d->DecayLength();
 608       Double_t cosp=d->CosPointingAngle();
 609       Double_t pt0=d->PtProng(0);
 610       Double_t pt1=d->PtProng(1);
 611       Double_t pt2=d->PtProng(2);
 612       Double_t sigvert=d->GetSigmaVert();
 613       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 614       Double_t dca=d->GetDCA();
 615       Double_t ptmax=0;
 616       for(Int_t i=0;i<3;i++){
 617         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 618       }
 619       fCosPHist[index]->Fill(cosp);
 620       fDLenHist[index]->Fill(dlen);
 621       fSigVertHist[index]->Fill(sigvert);
 622       fSumd02Hist[index]->Fill(sumD02);
 623       fPtMaxHist[index]->Fill(ptmax);
 624       fDCAHist[index]->Fill(dca);
 625       fPtProng0Hist[index]->Fill(pt0);
 626       fPtProng1Hist[index]->Fill(pt1);
 627       fPtProng2Hist[index]->Fill(pt2);
 628       if(isKKpi){
 629         Double_t massKK=d->InvMass2Prongs(0,1,321,321);
 630         Double_t massKp=d->InvMass2Prongs(1,2,321,211);
 631         fDalitz[index]->Fill(massKK,massKp);
 632         if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
 633         if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
 634         if(fReadMC && indexMCKKpi!=-1){
 635           fDalitz[indexMCKKpi]->Fill(massKK,massKp);
 636           if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
 637           if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
 638           fCosPHist[indexMCKKpi]->Fill(cosp);
 639           fDLenHist[indexMCKKpi]->Fill(dlen);
 640           fSigVertHist[indexMCKKpi]->Fill(sigvert);
 641           fSumd02Hist[indexMCKKpi]->Fill(sumD02);
 642           fPtMaxHist[indexMCKKpi]->Fill(ptmax);
 643           fPtCandHist[indexMCKKpi]->Fill(ptCand);
 644           fDCAHist[indexMCKKpi]->Fill(dca);
 645           fPtProng0Hist[indexMCKKpi]->Fill(pt0);
 646           fPtProng1Hist[indexMCKKpi]->Fill(pt1);
 647           fPtProng2Hist[indexMCKKpi]->Fill(pt2);
 648         }
 649       }
 650       if(ispiKK){
 651         Double_t massKK=d->InvMass2Prongs(1,2,321,321);
 652         Double_t massKp=d->InvMass2Prongs(0,1,211,321);
 653         fDalitz[index]->Fill(massKK,massKp);
 654         if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
 655         if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
 656
 657
 658         if(fReadMC && indexMCpiKK!=-1){
 659           fDalitz[indexMCpiKK]->Fill(massKK,massKp);
 660           if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
 661           if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
 662           fCosPHist[indexMCpiKK]->Fill(cosp);
 663           fDLenHist[indexMCpiKK]->Fill(dlen);
 664           fSigVertHist[indexMCpiKK]->Fill(sigvert);
 665           fSumd02Hist[indexMCpiKK]->Fill(sumD02);
 666           fPtMaxHist[indexMCpiKK]->Fill(ptmax);
 667           fPtCandHist[indexMCpiKK]->Fill(ptCand);
 668           fDCAHist[indexMCpiKK]->Fill(dca);
 669           fPtProng0Hist[indexMCpiKK]->Fill(pt0);
 670           fPtProng1Hist[indexMCpiKK]->Fill(pt1);
 671           fPtProng2Hist[indexMCpiKK]->Fill(pt2);
 672         }
 673       }
 674
 675
 676     }
 677
 678     Float_t tmp[26];
 679     if(fFillNtuple>0){
 680
 681       if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){
 682
 683         tmp[0]=Float_t(labDs);
 684         tmp[1]=Float_t(retCodeAnalysisCuts);
 685         tmp[2]=Float_t(pdgCode0);
 686         tmp[3]=d->PtProng(0);
 687         tmp[4]=d->PtProng(1);
 688         tmp[5]=d->PtProng(2);
 689         tmp[6]=d->Pt();
 690         tmp[7]=d->CosPointingAngle();
 691         tmp[8]=d->DecayLength();
 692         tmp[9]=d->Xv();
 693         tmp[10]=d->Yv();
 694         tmp[11]=d->Zv();
 695         tmp[12]=d->InvMassDsKKpi();
 696         tmp[13]=d->InvMassDspiKK();
 697         tmp[14]=d->GetSigmaVert();
 698         tmp[15]=d->Getd0Prong(0);
 699         tmp[16]=d->Getd0Prong(1);
 700         tmp[17]=d->Getd0Prong(2);
 701         tmp[18]=d->GetDCA();
 702         tmp[19]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 703         tmp[20]=d->InvMass2Prongs(0,1,321,321);
 704         tmp[21]=d->InvMass2Prongs(1,2,321,321);
 705         tmp[22]=d->CosPiDsLabFrameKKpi();
 706         tmp[23]=d->CosPiDsLabFramepiKK();
 707         tmp[24]=d->CosPiKPhiRFrameKKpi();
 708         tmp[25]=d->CosPiKPhiRFramepiKK();
 709
 710
 711         fNtupleDs->Fill(tmp);
 712         PostData(4,fNtupleDs);
 713       }
 714     }
 715
 716     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 717     if(recVtx)fProdCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
 718   }
 719
 720   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
 721   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
 722
 723   PostData(1,fOutput);
 724   PostData(3,fCounter);
 725
 726   return;
 727 }
 728
 729 //_________________________________________________________________
 730
 731 void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
 732 {
 733   // Terminate analysis
 734   //
 735   if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
 736   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 737   if (!fOutput) {
 738     printf("ERROR: fOutput not available\n");
 739     return;
 740   }
 741   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
 742   if(fHistNEvents){
 743     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
 744   }else{
 745     printf("ERROR: fHistNEvents not available\n");
 746     return;
 747   }
 748   return;
 749 }
 750