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