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