PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, 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 // Class AliAnalysisTaskSEDvsMultiplicity
  20 // AliAnalysisTaskSE for the D meson vs. multiplcity analysis
  21 // Authors: Renu Bala, Zaida Conesa del Valle, Francesco Prino
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TClonesArray.h>
  25 #include <TCanvas.h>
  26 #include <TList.h>
  27 #include <TString.h>
  28 #include <TDatabasePDG.h>
  29 #include <TH1F.h>
  30 #include <TH2F.h>
  31 #include <TH3F.h>
  32 #include <THnSparse.h>
  33 #include <TProfile.h>
  34 #include "AliAnalysisManager.h"
  35 #include "AliRDHFCuts.h"
  36 #include "AliRDHFCutsDplustoKpipi.h"
  37 #include "AliRDHFCutsDStartoKpipi.h"
  38 #include "AliRDHFCutsD0toKpi.h"
  39 #include "AliAODHandler.h"
  40 #include "AliAODEvent.h"
  41 #include "AliAODVertex.h"
  42 #include "AliAODTrack.h"
  43 #include "AliAODRecoDecayHF.h"
  44 #include "AliAODRecoCascadeHF.h"
  45 #include "AliAnalysisVertexingHF.h"
  46 #include "AliAnalysisTaskSE.h"
  47 #include "AliAnalysisTaskSEDvsMultiplicity.h"
  48 #include "AliNormalizationCounter.h"
  49 #include "AliVertexingHFUtils.h"
  50 ClassImp(AliAnalysisTaskSEDvsMultiplicity)
  51
  52
  53 //________________________________________________________________________
  54 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity():
  55 AliAnalysisTaskSE(),
  56   fOutput(0),
  57   fListCuts(0),
  58   fOutputCounters(0),
  59   fListProfiles(0),
  60   fHistNEvents(0),
  61   fHistNtrEta16vsNtrEta1(0),
  62   fHistNtrCorrEta1vsNtrRawEta1(0),
  63   fHistNtrVsZvtx(0),
  64   fHistNtrCorrVsZvtx(0),
  65   fHistNtrVsNchMC(0),
  66   fHistNtrCorrVsNchMC(0),
  67   fHistNtrVsNchMCPrimary(0),
  68   fHistNtrCorrVsNchMCPrimary(0),
  69   fHistNtrVsNchMCPhysicalPrimary(0),
  70   fHistNtrCorrVsNchMCPhysicalPrimary(0),
  71   fHistGenPrimaryParticlesInelGt0(0),
  72   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
  73   fHistNtrCorrEvSel(0),
  74   fHistNtrCorrEvWithCand(0),
  75   fHistNtrCorrEvWithD(0),
  76   fPtVsMassVsMult(0),
  77   fPtVsMassVsMultNoPid(0),
  78   fPtVsMassVsMultUncorr(0),
  79   fPtVsMassVsMultPart(0),
  80   fPtVsMassVsMultAntiPart(0),
  81   fUpmasslimit(1.965),
  82   fLowmasslimit(1.765),
  83   fNMassBins(200),
  84   fRDCutsAnalysis(0),
  85   fCounter(0),
  86   fCounterU(0),
  87   fDoImpPar(kFALSE),
  88   fNImpParBins(400),
  89   fLowerImpPar(-2000.),
  90   fHigherImpPar(2000.),
  91   fReadMC(kFALSE),
  92   fMCOption(0),
  93   fUseBit(kTRUE),
  94   fRefMult(9.5),
  95   fPdgMeson(411)
  96 {
  97    // Default constructor
  98   for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
  99   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
 100 }
 101
 102 //________________________________________________________________________
 103 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts):
 104   AliAnalysisTaskSE(name),
 105   fOutput(0),
 106   fListCuts(0),
 107   fOutputCounters(0),
 108   fListProfiles(0),
 109   fHistNEvents(0),
 110   fHistNtrEta16vsNtrEta1(0),
 111   fHistNtrCorrEta1vsNtrRawEta1(0),
 112   fHistNtrVsZvtx(0),
 113   fHistNtrCorrVsZvtx(0),
 114   fHistNtrVsNchMC(0),
 115   fHistNtrCorrVsNchMC(0),
 116   fHistNtrVsNchMCPrimary(0),
 117   fHistNtrCorrVsNchMCPrimary(0),
 118   fHistNtrVsNchMCPhysicalPrimary(0),
 119   fHistNtrCorrVsNchMCPhysicalPrimary(0),
 120   fHistGenPrimaryParticlesInelGt0(0),
 121   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
 122   fHistNtrCorrEvSel(0),
 123   fHistNtrCorrEvWithCand(0),
 124   fHistNtrCorrEvWithD(0),
 125   fPtVsMassVsMult(0),
 126   fPtVsMassVsMultNoPid(0),
 127   fPtVsMassVsMultUncorr(0),
 128   fPtVsMassVsMultPart(0),
 129   fPtVsMassVsMultAntiPart(0),
 130   fUpmasslimit(1.965),
 131   fLowmasslimit(1.765),
 132   fNMassBins(200),
 133   fRDCutsAnalysis(cuts),
 134   fCounter(0),
 135   fCounterU(0),
 136   fDoImpPar(kFALSE),
 137   fNImpParBins(400),
 138   fLowerImpPar(-2000.),
 139   fHigherImpPar(2000.),
 140   fReadMC(kFALSE),
 141   fMCOption(0),
 142   fUseBit(kTRUE),
 143   fRefMult(9.5),
 144   fPdgMeson(pdgMeson)
 145 {
 146   //
 147   // Standard constructor
 148   //
 149   for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
 150   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
 151   if(fPdgMeson==413){
 152     fNMassBins=200; // FIXME
 153     SetMassLimits(0.,0.2); // FIXME
 154   }else{
 155     fNMassBins=200;
 156     SetMassLimits(fPdgMeson,0.1);
 157   }
 158   // Default constructor
 159    // Output slot #1 writes into a TList container
 160   DefineOutput(1,TList::Class());  //My private output
 161   // Output slot #2 writes cut to private output
 162   DefineOutput(2,TList::Class());
 163   // Output slot #3 writes cut to private output
 164   DefineOutput(3,TList::Class());
 165   // Output slot #4 writes cut to private output
 166   DefineOutput(4,TList::Class());
 167 }
 168 //________________________________________________________________________
 169 AliAnalysisTaskSEDvsMultiplicity::~AliAnalysisTaskSEDvsMultiplicity()
 170 {
 171   //
 172   // Destructor
 173   //
 174   delete fOutput;
 175   delete fHistNEvents;
 176   delete fListCuts;
 177   delete fListProfiles;
 178   delete fRDCutsAnalysis;
 179   delete fCounter;
 180   delete fCounterU;
 181   for(Int_t i=0; i<4; i++) delete fMultEstimatorAvg[i];
 182   for(Int_t i=0; i<5; i++){
 183     delete fHistMassPtImpPar[i];
 184   }
 185 }
 186
 187 //_________________________________________________________________
 188 void  AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Double_t lowlimit, Double_t uplimit){
 189   // set invariant mass limits
 190   if(uplimit>lowlimit){
 191     fLowmasslimit = lowlimit;
 192     fUpmasslimit = uplimit;
 193   }else{
 194     AliError("Wrong mass limits: upper value should be larger than lower one");
 195   }
 196 }
 197 //_________________________________________________________________
 198 void  AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Int_t pdg, Double_t range){
 199   // set invariant mass limits
 200   Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
 201   SetMassLimits(mass-range,mass+range);
 202 }
 203 //________________________________________________________________________
 204 void AliAnalysisTaskSEDvsMultiplicity::Init(){
 205   //
 206   // Initialization
 207   //
 208   printf("AnalysisTaskSEDvsMultiplicity::Init() \n");
 209
 210   fListCuts=new TList();
 211
 212   if(fPdgMeson==411){
 213      AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
 214     copycut->SetName("AnalysisCutsDplus");
 215     fListCuts->Add(copycut);
 216   }else if(fPdgMeson==421){
 217     AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
 218     copycut->SetName("AnalysisCutsDzero");
 219     fListCuts->Add(copycut);
 220   }else if(fPdgMeson==413){
 221      AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
 222     copycut->SetName("AnalysisCutsDStar");
 223     fListCuts->Add(copycut);
 224   }
 225   PostData(2,fListCuts);
 226
 227   fListProfiles = new TList();
 228   fListProfiles->SetOwner();
 229   TString period[4]={"LHC10b","LHC10c","LHC10d","LHC10e"};
 230   for(Int_t i=0; i<4; i++){
 231     if(fMultEstimatorAvg[i]){
 232       TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
 233       hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
 234       fListProfiles->Add(hprof);
 235     }
 236   }
 237   PostData(4,fListProfiles);
 238
 239   return;
 240 }
 241
 242 //________________________________________________________________________
 243 void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
 244 {
 245   // Create the output container
 246   //
 247   if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
 248
 249   // Several histograms are more conveniently managed in a TList
 250   fOutput = new TList();
 251   fOutput->SetOwner();
 252   fOutput->SetName("OutputHistos");
 253
 254   fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel","Tracklets multiplicity for selected events; Tracklets ; Entries",200,0.,200.);
 255   fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",200,0.,200.);// Total multiplicity
 256   fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",200,0.,200.); //
 257   fHistNtrEta16vsNtrEta1 = new TH2F("hNtrEta16vsNtrEta1","Uncorrected Eta1.6 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<1.6",200,-0.5,199.5,200,-0.5,199.5); //eta 1.6 vs eta 1.0 histogram
 258   fHistNtrCorrEta1vsNtrRawEta1 = new TH2F("hNtrCorrEta1vsNtrRawEta1","Corrected Eta1 vs Eta1.0; Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",200,-0.,200.,200,-0.5,199.5); //eta 1.6 vs eta 1.0 histogram
 259   fHistNtrVsZvtx = new TH2F("hNtrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,0,200.); //
 260   fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,0,200.); //
 261
 262   fHistNtrVsNchMC = new TH2F("hNtrVsNchMC","Ntracklet vs NchMC; Nch;N_{tracklet};",200,0,200,200,0,200.); //
 263   fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC","Ntracklet vs Nch; Nch;N_{tracklet};",200,0,200,200,0,200.); //
 264
 265   fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch (Primary);N_{tracklet};",200,0,200,200,0,200.); //
 266   fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch(Primary) ;N_{tracklet};",200,0,200,200,0,200.); //
 267
 268   fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",200,0,200,200,0,200.); //
 269   fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",200,0,200,200,0,200.); //
 270
 271   fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",200,-0.5,199.5);
 272
 273   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary = new TH3F("fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary", "MC: Nch (Physical Primary) vs Nch (Primary) vs Nch (Generated); Nch (Generated); Nch (Primary); Nch (Physical Primary)",200,0.,200.,200,0.,200.,200,0.,200.);
 274
 275   fHistNtrCorrEvSel->Sumw2();
 276   fHistNtrCorrEvWithCand->Sumw2();
 277   fHistNtrCorrEvWithD->Sumw2();
 278   fHistGenPrimaryParticlesInelGt0->Sumw2();
 279   fOutput->Add(fHistNtrCorrEvSel);
 280   fOutput->Add(fHistNtrCorrEvWithCand);
 281   fOutput->Add(fHistNtrCorrEvWithD);
 282   fOutput->Add(fHistNtrEta16vsNtrEta1);
 283   fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1);
 284   fOutput->Add(fHistNtrVsZvtx);
 285   fOutput->Add(fHistNtrCorrVsZvtx);
 286
 287   fOutput->Add(fHistNtrVsNchMC);
 288   fOutput->Add(fHistNtrCorrVsNchMC);
 289   fOutput->Add(fHistNtrVsNchMCPrimary);
 290   fOutput->Add(fHistNtrCorrVsNchMCPrimary);
 291   fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
 292   fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
 293   fOutput->Add(fHistGenPrimaryParticlesInelGt0);
 294   fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
 295
 296
 297   fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
 298   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
 299   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
 300   fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
 301   fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
 302   fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
 303   fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
 304   fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
 305   fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
 306   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
 307   fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
 308   fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
 309   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 310   fHistNEvents->Sumw2();
 311   fHistNEvents->SetMinimum(0);
 312   fOutput->Add(fHistNEvents);
 313
 314   fPtVsMassVsMult=new TH3F("hPtVsMassvsMult", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 315
 316   fPtVsMassVsMultNoPid=new TH3F("hPtVsMassvsMultNoPid", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 317
 318   fPtVsMassVsMultUncorr=new TH3F("hPtVsMassvsMultUncorr", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 319
 320   fPtVsMassVsMultPart=new TH3F("hPtVsMassvsMultPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 321
 322   fPtVsMassVsMultAntiPart=new TH3F("hPtVsMassvsMultAntiPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 323
 324   fOutput->Add(fPtVsMassVsMult);
 325   fOutput->Add(fPtVsMassVsMultUncorr);
 326   fOutput->Add(fPtVsMassVsMultNoPid);
 327   fOutput->Add(fPtVsMassVsMultPart);
 328   fOutput->Add(fPtVsMassVsMultAntiPart);
 329
 330   if(fDoImpPar) CreateImpactParameterHistos();
 331
 332   fCounter = new AliNormalizationCounter("NormCounterCorrMult");
 333   fCounter->SetStudyMultiplicity(kTRUE,1.);
 334   fCounter->Init();
 335
 336   fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
 337   fCounterU->SetStudyMultiplicity(kTRUE,1.);
 338   fCounterU->Init();
 339
 340   fOutputCounters = new TList();
 341   fOutputCounters->SetOwner();
 342   fOutputCounters->SetName("OutputCounters");
 343   fOutputCounters->Add(fCounter);
 344   fOutputCounters->Add(fCounterU);
 345
 346   PostData(1,fOutput);
 347   PostData(2,fListCuts);
 348   PostData(3,fOutputCounters);
 349   PostData(4,fListProfiles);
 350
 351
 352
 353   return;
 354 }
 355
 356
 357
 358 //________________________________________________________________________
 359 void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
 360 {
 361   // Execute analysis for current event:
 362   // heavy flavor candidates association to MC truth
 363
 364   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 365
 366   //  AliAODTracklets* tracklets = aod->GetTracklets();
 367   //Int_t ntracklets = tracklets->GetNumberOfTracklets();
 368
 369
 370   TClonesArray *arrayCand = 0;
 371   TString arrayName="";
 372   UInt_t pdgDau[3];
 373   Int_t nDau=0;
 374   Int_t selbit=0;
 375   if(fPdgMeson==411){
 376     arrayName="Charm3Prong";
 377     pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
 378     nDau=3;
 379     selbit=AliRDHFCuts::kDplusCuts;
 380   }else if(fPdgMeson==421){
 381     arrayName="D0toKpi";
 382     pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
 383     nDau=2;
 384     selbit=AliRDHFCuts::kD0toKpiCuts;
 385   }else if(fPdgMeson==413){
 386     arrayName="Dstar";
 387     pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=211;
 388     nDau=3;
 389     selbit=AliRDHFCuts::kDstarCuts;
 390   }
 391
 392   if(!aod && AODEvent() && IsStandardAOD()) {
 393     // In case there is an AOD handler writing a standard AOD, use the AOD
 394     // event in memory rather than the input (ESD) event.
 395     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 396      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 397      // have to taken from the AOD event hold by the AliAODExtension
 398     AliAODHandler* aodHandler = (AliAODHandler*)
 399       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 400     if(aodHandler->GetExtensions()) {
 401       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 402       AliAODEvent *aodFromExt = ext->GetAOD();
 403       arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
 404     }
 405   } else if(aod) {
 406     arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
 407   }
 408
 409   if(!aod || !arrayCand) {
 410     printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
 411     return;
 412   }
 413
 414
 415
 416   // fix for temporary bug in ESDfilter
 417   // the AODs with null vertex pointer didn't pass the PhysSel
 418   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
 419
 420   Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
 421   Int_t countTr=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
 422
 423   fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countTreta1);
 424   fHistNEvents->Fill(0); // count event
 425
 426   Double_t countTreta1corr=countTreta1;
 427   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 428   if(vtx1){
 429     if(vtx1->GetNContributors()>0){
 430       fHistNEvents->Fill(1);
 431       TProfile* estimatorAvg = GetEstimatorHistogram(aod);
 432       if(estimatorAvg){
 433         countTreta1corr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult);
 434       }
 435     }
 436   }
 437
 438   fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countTreta1corr);
 439
 440   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
 441
 442   if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
 443   if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
 444   if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
 445   if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
 446
 447
 448   if(!isEvSel)return;
 449   fHistNtrEta16vsNtrEta1->Fill(countTreta1,countTr);
 450   fHistNtrCorrEta1vsNtrRawEta1->Fill(countTreta1,countTreta1corr);
 451   if(vtx1){
 452     fHistNtrVsZvtx->Fill(vtx1->GetZ(),countTreta1);
 453     fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countTreta1corr);
 454   }
 455
 456   TClonesArray *arrayMC=0;
 457   AliAODMCHeader *mcHeader=0;
 458
 459   // load MC particles
 460   if(fReadMC){
 461
 462     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 463     if(!arrayMC) {
 464       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
 465       return;
 466     }
 467     // load MC header
 468     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 469     if(!mcHeader) {
 470       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
 471       return;
 472      }
 473
 474
 475     Int_t nChargedMC=AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.0,1.0);
 476     Int_t nChargedMCPrimary=AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.0,1.0);
 477     Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.0,1.0);
 478     if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
 479       fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary);
 480     }
 481     fHistNtrVsNchMC->Fill(nChargedMC,countTreta1);
 482     fHistNtrCorrVsNchMC->Fill(nChargedMC,countTreta1corr);
 483
 484     fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1);
 485     fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1corr);
 486
 487     fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1);
 488     fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1corr);
 489
 490     fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary);
 491   }
 492
 493   Int_t nCand = arrayCand->GetEntriesFast();
 494   Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
 495   Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 496   Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
 497   Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 498
 499   for (Int_t iCand = 0; iCand < nCand; iCand++) {
 500     AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
 501     fHistNEvents->Fill(7);
 502     if(fUseBit && !d->HasSelectionBit(selbit)){
 503       fHistNEvents->Fill(8);
 504       continue;
 505     }
 506
 507     Double_t ptCand = d->Pt();
 508     Double_t rapid=d->Y(fPdgMeson);
 509     Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
 510     if(!isFidAcc) continue;
 511     Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
 512     Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
 513     if(passTopolCuts==0) continue;
 514     nSelectedNoPID++;
 515     fHistNEvents->Fill(9);
 516     if(passAllCuts){
 517       nSelectedPID++;
 518       fHistNEvents->Fill(10);
 519     }
 520     Bool_t isPrimary=kTRUE;
 521     Int_t labD=-1;
 522     Double_t trueImpParXY=9999.;
 523     Double_t impparXY=d->ImpParXY()*10000.;
 524     Double_t dlen=0.1; //FIXME
 525     Double_t mass[2];
 526     if(fPdgMeson==411){
 527       mass[0]=d->InvMass(nDau,pdgDau);
 528       mass[1]=-1.;
 529       if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
 530     }else if(fPdgMeson==421){
 531       UInt_t pdgdaughtersD0[2]={211,321};//pi,K
 532       UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
 533       mass[0]=d->InvMass(2,pdgdaughtersD0);
 534       mass[1]=d->InvMass(2,pdgdaughtersD0bar);
 535       if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
 536     }else if(fPdgMeson==413){
 537       // FIXME
 538       AliAODRecoCascadeHF* temp = (AliAODRecoCascadeHF*)d;
 539       mass[0]=temp->DeltaInvMass();
 540       mass[1]=-1.;
 541       if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.001) nSelectedInMassPeak++; //1 MeV for now... FIXME
 542     }
 543     for(Int_t iHyp=0; iHyp<2; iHyp++){
 544       if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
 545       Double_t invMass=mass[iHyp];
 546       Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,countTreta1corr};
 547
 548       if(fReadMC){
 549
 550         labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
 551         Bool_t fillHisto=fDoImpPar;
 552         if(labD>=0){
 553           AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
 554           Int_t code=partD->GetPdgCode();
 555           if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
 556           if(code<0 && iHyp==0) fillHisto=kFALSE;
 557           if(code>0 && iHyp==1) fillHisto=kFALSE;
 558           if(!isPrimary){
 559             if(fPdgMeson==411){
 560               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
 561             }else if(fPdgMeson==421){
 562               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
 563             }else if(fPdgMeson==413){
 564               trueImpParXY=0.; /// FIXME
 565             }
 566             Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,countTreta1corr};
 567             if(fillHisto && passAllCuts){
 568               fHistMassPtImpPar[2]->Fill(arrayForSparse);
 569               fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
 570             }
 571           }else{
 572             if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
 573           }
 574         }else{
 575           if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
 576         }
 577         if(fPdgMeson==421){
 578           if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
 579           if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
 580         }
 581       }
 582
 583       if(fPdgMeson==421){
 584         if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
 585         if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
 586       }
 587
 588       fPtVsMassVsMultNoPid->Fill(countTreta1corr,invMass,ptCand);
 589       if(passAllCuts){
 590         fPtVsMassVsMult->Fill(countTreta1corr,invMass,ptCand);
 591         fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand);
 592         // Add separation between part antipart
 593         if(fPdgMeson==411){
 594           if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
 595           else fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
 596         }else if(fPdgMeson==421){
 597           if(passTopolCuts&1) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
 598           if(passTopolCuts&2) fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
 599         }else if(fPdgMeson==413){
 600           // FIXME ADD Dstar!!!!!!!!
 601         }
 602
 603         if(fDoImpPar){
 604           fHistMassPtImpPar[0]->Fill(arrayForSparse);
 605         }
 606
 607       }
 608
 609     }
 610   }
 611   fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
 612   fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
 613   fHistNtrCorrEvSel->Fill(countTreta1corr);
 614   if(nSelectedPID>0) fHistNtrCorrEvWithCand->Fill(countTreta1corr);
 615   if(nSelectedInMassPeak>0) fHistNtrCorrEvWithD->Fill(countTreta1corr);
 616
 617   PostData(1,fOutput);
 618   PostData(2,fListCuts);
 619   PostData(3,fOutputCounters);
 620
 621   return;
 622 }
 623 //________________________________________________________________________
 624 void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
 625   // Histos for impact paramter study
 626   // mass . pt , impact parameter , decay length , multiplicity
 627
 628   Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
 629   Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
 630   Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
 631
 632   fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
 633                                         "Mass vs. pt vs.imppar - All",
 634                                         5,nbins,xmin,xmax);
 635   fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
 636                                         "Mass vs. pt vs.imppar - promptD",
 637                                         5,nbins,xmin,xmax);
 638   fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
 639                                         "Mass vs. pt vs.imppar - DfromB",
 640                                         5,nbins,xmin,xmax);
 641   fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
 642                                         "Mass vs. pt vs.true imppar -DfromB",
 643                                         5,nbins,xmin,xmax);
 644   fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
 645                                         "Mass vs. pt vs.imppar - backgr.",
 646                                         5,nbins,xmin,xmax);
 647   for(Int_t i=0; i<5;i++){
 648     fOutput->Add(fHistMassPtImpPar[i]);
 649   }
 650 }
 651
 652 //________________________________________________________________________
 653 void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
 654 {
 655   // Terminate analysis
 656   //
 657   if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
 658
 659   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 660   if (!fOutput) {
 661     printf("ERROR: fOutput not available\n");
 662     return;
 663   }
 664
 665   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 666   if(!fHistNEvents){
 667     printf("ERROR: fHistNEvents not available\n");
 668     return;
 669   }
 670   printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
 671
 672   return;
 673 }
 674 //_________________________________________________________________________________________________
 675 Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
 676   //
 677   // checking whether the mother of the particles come from a charm or a bottom quark
 678   //
 679
 680   Int_t pdgGranma = 0;
 681   Int_t mother = 0;
 682   mother = mcPartCandidate->GetMother();
 683   Int_t istep = 0;
 684   Int_t abspdgGranma =0;
 685   Bool_t isFromB=kFALSE;
 686   Bool_t isQuarkFound=kFALSE;
 687   while (mother >0 ){
 688     istep++;
 689     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
 690     if (mcGranma){
 691       pdgGranma = mcGranma->GetPdgCode();
 692       abspdgGranma = TMath::Abs(pdgGranma);
 693       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 694         isFromB=kTRUE;
 695       }
 696       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 697       mother = mcGranma->GetMother();
 698     }else{
 699       AliError("Failed casting the mother particle!");
 700       break;
 701     }
 702   }
 703
 704   if(isFromB) return 5;
 705   else return 4;
 706 }
 707
 708
 709
 710 //____________________________________________________________________________
 711 TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
 712   // Get Estimator Histogram from period event->GetRunNumber();
 713   //
 714   // If you select SPD tracklets in |eta|<1 you should use type == 1
 715   //
 716
 717   Int_t runNo  = event->GetRunNumber();
 718   Int_t period = -1;   // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
 719   if(runNo>114930 && runNo<117223) period = 0;
 720   if(runNo>119158 && runNo<120830) period = 1;
 721   if(runNo>122373 && runNo<126438) period = 2;
 722   if(runNo>127711 && runNo<130841) period = 3;
 723   if(period<0 || period>3) return 0;
 724
 725   return fMultEstimatorAvg[period];
 726 }