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