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