PWG3/vertexingHF/AliAnalysisTaskSEHFQA.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-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 /////////////////////////////////////////////////////////////
  17 //
  18 // AliAnalysisTaskSE for HF quality assurance
  19 //
  20 // Author: Chiara Bianchin, chiara.bianchin@pd.infn.it
  21 /////////////////////////////////////////////////////////////
  22
  23 #include <Riostream.h>
  24 #include <TClonesArray.h>
  25 #include <TCanvas.h>
  26 #include <TList.h>
  27 #include <TH1F.h>
  28 #include <TH2F.h>
  29 #include <TDatabasePDG.h>
  30
  31 #include <AliAnalysisDataSlot.h>
  32 #include <AliAnalysisDataContainer.h>
  33 #include "AliAnalysisManager.h"
  34 #include "AliESDtrack.h"
  35 #include "AliVertexerTracks.h"
  36 #include "AliPID.h"
  37 #include "AliTPCPIDResponse.h"
  38 #include "AliAODHandler.h"
  39 #include "AliAODEvent.h"
  40 #include "AliAODVertex.h"
  41 #include "AliAODTrack.h"
  42 #include "AliAODMCParticle.h"
  43 #include "AliAODMCHeader.h"
  44 #include "AliAODRecoDecayHF2Prong.h"
  45 #include "AliAODRecoCascadeHF.h"
  46 #include "AliAnalysisVertexingHF.h"
  47 #include "AliAnalysisTaskSE.h"
  48 #include "AliCounterCollection.h"
  49 #include "AliRDHFCuts.h"
  50 #include "AliRDHFCutsDplustoKpipi.h"
  51 #include "AliRDHFCutsD0toKpipipi.h"
  52 #include "AliRDHFCutsDstoKKpi.h"
  53 #include "AliRDHFCutsDStartoKpipi.h"
  54 #include "AliRDHFCutsD0toKpi.h"
  55 #include "AliRDHFCutsLctopKpi.h"
  56
  57 #include "AliAnalysisTaskSEHFQA.h"
  58
  59
  60 ClassImp(AliAnalysisTaskSEHFQA)
  61
  62 //____________________________________________________________________________
  63
  64 AliAnalysisTaskSEHFQA::AliAnalysisTaskSEHFQA():AliAnalysisTaskSE(),
  65   fNEntries(0x0),
  66   fOutputPID(0x0),
  67   fOutputTrack(0x0),
  68   fOutputCounters(0x0),
  69   fOutputCheckCentrality(0x0),
  70   fDecayChannel(AliAnalysisTaskSEHFQA::kD0toKpi),
  71   fCuts(0x0),
  72   fEstimator(AliRDHFCuts::kCentTRK),
  73   fReadMC(kFALSE),
  74   fSimpleMode(kFALSE)
  75 {
  76   //default constructor
  77 }
  78
  79 //____________________________________________________________________________
  80 AliAnalysisTaskSEHFQA::AliAnalysisTaskSEHFQA(const char *name, AliAnalysisTaskSEHFQA::DecChannel ch,AliRDHFCuts* cuts):
  81   AliAnalysisTaskSE(name),
  82   fNEntries(0x0),
  83   fOutputPID(0x0),
  84   fOutputTrack(0x0),
  85   fOutputCounters(0x0),
  86   fOutputCheckCentrality(0x0),
  87   fDecayChannel(ch),
  88   fCuts(0x0),
  89   fEstimator(AliRDHFCuts::kCentTRK),
  90   fReadMC(kFALSE),
  91   fSimpleMode(kFALSE)
  92 {
  93   //constructor
  94
  95   //SetCutObject(cuts);
  96   fCuts=cuts;
  97
  98   // Output slot #1 writes into a TH1F container (number of events)
  99   DefineOutput(1,TH1F::Class());  //My private output
 100   // Output slot #2 writes into a TList container (PID)
 101   DefineOutput(2,TList::Class());  //My private output
 102   // Output slot #3 writes into a TList container (Tracks)
 103   DefineOutput(3,TList::Class());  //My private output
 104   // Output slot #4 writes into a AliRDHFCuts container (cuts)
 105   switch(fDecayChannel){
 106   case 0:
 107     DefineOutput(4,AliRDHFCutsDplustoKpipi::Class());  //My private output
 108     break;
 109   case 1:
 110     DefineOutput(4,AliRDHFCutsD0toKpi::Class());  //My private output
 111     break;
 112   case 2:
 113     DefineOutput(4,AliRDHFCutsDStartoKpipi::Class());  //My private output
 114     break;
 115   case 3:
 116     DefineOutput(4,AliRDHFCutsDstoKKpi::Class());  //My private output
 117     break;
 118   case 4:
 119     DefineOutput(4,AliRDHFCutsD0toKpipipi::Class());  //My private output
 120     break;
 121   case 5:
 122     DefineOutput(4,AliRDHFCutsLctopKpi::Class());  //My private output
 123     break;
 124   }
 125   // Output slot #5 writes into a TList container (AliCounterCollection)
 126   DefineOutput(5,TList::Class());  //My private output
 127   // Output slot #6 writes into a TList container (TH1F)
 128   DefineOutput(6,TList::Class());  //My private output
 129 }
 130
 131 //___________________________________________________________________________
 132 AliAnalysisTaskSEHFQA::~AliAnalysisTaskSEHFQA()
 133 {
 134   //destructor
 135   delete fNEntries;
 136
 137   delete fOutputPID;
 138
 139   delete fOutputTrack;
 140
 141   delete fOutputCounters;
 142
 143   delete fOutputCheckCentrality;
 144
 145 }
 146
 147 //___________________________________________________________________________
 148 void AliAnalysisTaskSEHFQA::Init(){
 149
 150   //initialization
 151   if(fDebug > 1) printf("AnalysisTaskSEHFQA::Init() \n");
 152
 153   switch(fDecayChannel){
 154   case 0:
 155     {
 156       AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
 157       // Post the data
 158       PostData(4,copycut);
 159     }
 160     break;
 161   case 1:
 162     {
 163       AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 164       // Post the data
 165       PostData(4,copycut);
 166     }
 167     break;
 168   case 2:
 169     {
 170       AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 171       // Post the data
 172       PostData(4,copycut);
 173     }
 174     break;
 175   case 3:
 176     {
 177       AliRDHFCutsDstoKKpi* copycut=new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
 178       // Post the data
 179       PostData(4,copycut);
 180     }
 181     break;
 182   case 4:
 183     {
 184       AliRDHFCutsD0toKpipipi* copycut=new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
 185       // Post the data
 186       PostData(4,copycut);
 187     }
 188     break;
 189   case 5:
 190     {
 191       AliRDHFCutsLctopKpi* copycut=new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
 192       // Post the data
 193       PostData(4,copycut);
 194     }
 195     break;
 196
 197   default:
 198     return;
 199   }
 200
 201
 202
 203 }
 204
 205 //___________________________________________________________________________
 206 void AliAnalysisTaskSEHFQA::UserCreateOutputObjects()
 207 {
 208
 209   //create the output container
 210   if(fDebug > 1) printf("AnalysisTaskSEHFQA::UserCreateOutputObjects() \n");
 211
 212   //count events
 213
 214   fNEntries=new TH1F(GetOutputSlot(1)->GetContainer()->GetName(), "Counts the number of events", 9,-0.5,8.5);
 215   fNEntries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 216   fNEntries->GetXaxis()->SetBinLabel(2,"Pile-up Rej");
 217   fNEntries->GetXaxis()->SetBinLabel(3,"No VertexingHF");
 218   fNEntries->GetXaxis()->SetBinLabel(4,"nCandidates(AnCuts)");
 219   fNEntries->GetXaxis()->SetBinLabel(5,"EventsWithGoodVtx");
 220   fNEntries->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
 221   if(fReadMC){
 222     fNEntries->GetXaxis()->SetBinLabel(7,"MC Cand from c");
 223     fNEntries->GetXaxis()->SetBinLabel(8,"MC Cand from b");
 224     fNEntries->GetXaxis()->SetBinLabel(9,"N fakes Trks");
 225   } else{
 226     fNEntries->GetXaxis()->SetBinLabel(7,"N candidates");
 227   }
 228
 229   fNEntries->GetXaxis()->SetNdivisions(1,kFALSE);
 230
 231   //PID
 232
 233   fOutputPID=new TList();
 234   fOutputPID->SetOwner();
 235   fOutputPID->SetName(GetOutputSlot(2)->GetContainer()->GetName());
 236
 237   //TOF pid
 238   TString hname="hTOFsig";
 239   TH1F* hTOFsig=new TH1F(hname.Data(),"Distribution of TOF signal;TOF time [ps];Entries", 100, -2.e3,40.e3);
 240
 241   hname="hTOFtime";
 242   TH1F* hTOFtime=new TH1F(hname.Data(),"Distribution of TOF time Kaon;TOF time(Kaon) [ps];Entries", 1000, 0.,50000.);
 243
 244   hname="hTOFtimeKaonHyptime";
 245   TH2F* hTOFtimeKaonHyptime=new TH2F(hname.Data(),"TOFtime - timeHypothesisForKaon;p[GeV/c];TOFtime - timeHypothesisForKaon [ps]",200,0.,4.,1000,-20000.,20000.);
 246
 247   hname="hTOFtimeKaonHyptimeAC";
 248   TH2F* hTOFtimeKaonHyptimeAC=new TH2F(hname.Data(),"TOFtime - timeHypothesisForKaon;p[GeV/c];TOFtime - timeHypothesisForKaon [ps]",200,0.,4.,1000,-20000.,20000.);
 249
 250   hname="hTOFsigmaK160";
 251   TH2F* hTOFsigmaK160=new TH2F(hname.Data(),"(TOFsignal-timeK)/160 ps;p[GeV/c];(TOFsignal-timeK)/160 ps",200,0.,4.,100,-5,5);
 252
 253   hname="hTOFsigmaPion160";
 254   TH2F* hTOFsigmaPion160=new TH2F(hname.Data(),"(TOFsignal-time#pi)/160 ps;p[GeV/c];(TOFsignal-time#pi)/160 ps",200,0.,4.,100,-5,5);
 255
 256   hname="hTOFsigmaProton160";
 257   TH2F* hTOFsigmaProton160=new TH2F(hname.Data(),"(TOFsignal-timep)/160 ps;p[GeV/c];(TOFsignal-time p)/160 ps",200,0.,4.,100,-5,5);
 258
 259   hname="hTOFsigmaKSigPid";
 260   TH2F* hTOFsigmaKSigPid=new TH2F(hname.Data(),"(TOFsignal-timeK)/tofSigPid;p[GeV/c];(TOFsignal-timeK)/tofSigPid",200,0.,4.,100,-5,5);
 261
 262   hname="hTOFsigmaPionSigPid";
 263   TH2F* hTOFsigmaPionSigPid=new TH2F(hname.Data(),"(TOFsignal-time#pi)/tofSigPid;p[GeV/c];(TOFsignal-time#pi)/tofSigPid",200,0.,4.,100,-5,5);
 264
 265   hname="hTOFsigmaProtonSigPid";
 266   TH2F* hTOFsigmaProtonSigPid=new TH2F(hname.Data(),"(TOFsignal-timep)/tofSigPid;p[GeV/c];(TOFsignal-time p)/tofSigPid",200,0.,4.,100,-5,5);
 267
 268   hname="hTOFsigPid3sigPion";
 269   TH1F* hTOFsigPid3sigPion=new TH1F(hname.Data(),"TOF PID resolution (#pi) [ps]",500,0.,1000.);
 270
 271   hname="hTOFsigPid3sigKaon";
 272   TH1F* hTOFsigPid3sigKaon=new TH1F(hname.Data(),"TOF PID resolution (K) [ps]",500,0.,1000.);
 273
 274   hname="hTOFsigPid3sigProton";
 275   TH1F* hTOFsigPid3sigProton=new TH1F(hname.Data(),"TOF PID resolution (p) [ps]",500,0.,1000.);
 276
 277
 278   //TPC pid
 279   hname="hTPCsig";
 280   TH1F* hTPCsig=new TH1F(hname.Data(),"Distribution of TPC signal;TPC sig;Entries", 100, 35.,100.);
 281
 282   hname="hTPCsigvsp";
 283   TH2F* hTPCsigvsp=new TH2F(hname.Data(),"TPCsig vs p;TPC p[GeV/c];TPCsig",200,0.,4.,1000,35.,100.);
 284
 285   hname="hTPCsigvspAC";
 286   TH2F* hTPCsigvspAC=new TH2F(hname.Data(),"TPCsig vs p;TPCp[GeV/c];TPCsig",200,0.,4.,1000,35.,100.);
 287
 288   hname="hTPCsigmaK";
 289   TH2F* hTPCsigmaK=new TH2F(hname.Data(),"TPC Sigma for K as a function of momentum;p[GeV/c];Sigma Kaon",200,0.,4.,200,-5,5);
 290
 291   hname="hTPCsigmaPion";
 292   TH2F* hTPCsigmaPion=new TH2F(hname.Data(),"TPC Sigma for #pi as a function of momentum;p[GeV/c];Sigma #pi",200,0.,4.,200,-5,5);
 293
 294   hname="hTPCsigmaProton";
 295   TH2F* hTPCsigmaProton=new TH2F(hname.Data(),"TPC Sigma for proton as a function of momentum;p[GeV/c];Sigma Proton",200,0.,4.,200,-5,5);
 296
 297   fOutputPID->Add(hTOFsig);
 298   fOutputPID->Add(hTPCsig);
 299   fOutputPID->Add(hTOFtime);
 300   fOutputPID->Add(hTOFtimeKaonHyptime);
 301   fOutputPID->Add(hTOFtimeKaonHyptimeAC);
 302   fOutputPID->Add(hTOFsigmaK160);
 303   fOutputPID->Add(hTOFsigmaPion160);
 304   fOutputPID->Add(hTOFsigmaProton160);
 305   fOutputPID->Add(hTOFsigmaKSigPid);
 306   fOutputPID->Add(hTOFsigmaPionSigPid);
 307   fOutputPID->Add(hTOFsigmaProtonSigPid);
 308   fOutputPID->Add(hTOFsigPid3sigPion);
 309   fOutputPID->Add(hTOFsigPid3sigKaon);
 310   fOutputPID->Add(hTOFsigPid3sigProton);
 311   fOutputPID->Add(hTPCsigvsp);
 312   fOutputPID->Add(hTPCsigvspAC);
 313   fOutputPID->Add(hTPCsigmaK);
 314   fOutputPID->Add(hTPCsigmaPion);
 315   fOutputPID->Add(hTPCsigmaProton);
 316
 317   //quality of the tracks
 318
 319   fOutputTrack=new TList();
 320   fOutputTrack->SetOwner();
 321   fOutputTrack->SetName(GetOutputSlot(3)->GetContainer()->GetName());
 322
 323   hname="hnClsITS";
 324   TH1F* hnClsITS=new TH1F(hname.Data(),"Distribution of number of ITS clusters;nITScls;Entries",7,-0.5,6.5);
 325
 326   hname="hnClsITS-SA";
 327   TH1F* hnClsITSSA=new TH1F(hname.Data(),"Distribution of number of ITS clusters(ITS-SA);nITScls;Entries",7,-0.5,6.5);
 328
 329   hname="hnClsSPD";
 330   TH1F* hnClsSPD=new TH1F(hname.Data(),"Distribution of number of SPD clusters;nSPDcls;Entries",3,-0.5,2.5);
 331
 332   hname="hptGoodTr";
 333   TH1F* hptGoodTr=new TH1F(hname.Data(),"Pt distribution of 'good' tracks;p_{t}[GeV];Entries/0.05 GeV/c",400,0.,20.);
 334   hptGoodTr->SetTitleOffset(1.3,"Y");
 335
 336   hname="hdistrGoodTr";
 337   TH1F* hdistrGoodTr=new TH1F(hname.Data(),"Distribution of number of 'good' tracks per event;no.good-tracks/ev;Entries",4000,-0.5,3999.5);
 338   hdistrGoodTr->SetTitleOffset(1.3,"Y");
 339
 340   hname="hd0";
 341   TH1F* hd0=new TH1F(hname.Data(),"Impact parameter distribution of 'good' tracks;d_{0}[cm];Entries/10^{3} cm",200,-0.1,0.1);
 342
 343   fOutputTrack->Add(hnClsITS);
 344   fOutputTrack->Add(hnClsITSSA);
 345   fOutputTrack->Add(hnClsSPD);
 346   fOutputTrack->Add(hptGoodTr);
 347   fOutputTrack->Add(hdistrGoodTr);
 348   fOutputTrack->Add(hd0);
 349
 350   if(fCuts->GetUseCentrality()){
 351
 352     //Centrality (Counters)
 353     fOutputCounters=new TList();
 354     fOutputCounters->SetOwner();
 355     fOutputCounters->SetName(GetOutputSlot(5)->GetContainer()->GetName());
 356
 357     AliCounterCollection *stdEstimator=new AliCounterCollection("stdEstimator");
 358     stdEstimator->AddRubric("run",500000);
 359     stdEstimator->AddRubric("centralityclass","0_10/10_20/20_30/30_40/40_50/50_60/60_70/70_80/80_90/90_100");
 360     stdEstimator->Init();
 361     AliCounterCollection *secondEstimator=new AliCounterCollection("secondEstimator");
 362     secondEstimator->AddRubric("run",500000);
 363     secondEstimator->AddRubric("centralityclass","0_10/10_20/20_30/30_40/40_50/50_60/60_70/70_80/80_90/90_100/-990_-980");
 364     secondEstimator->Init();
 365
 366     fOutputCounters->Add(stdEstimator);
 367     fOutputCounters->Add(secondEstimator);
 368
 369     //Centrality (Checks)
 370     fOutputCheckCentrality=new TList();
 371     fOutputCheckCentrality->SetOwner();
 372     fOutputCheckCentrality->SetName(GetOutputSlot(6)->GetContainer()->GetName());
 373
 374     hname="hNtrackletsIn";
 375     TH1F* hNtrackletsIn=new TH1F(hname.Data(),"Number of tracklets in Centrality range;ntracklets;Entries",5000,-0.5,4999.5);
 376
 377     hname="hMultIn";
 378     TH1F* hMultIn=new TH1F(hname.Data(),"Multiplicity;multiplicity in Centrality range;Entries",10000,-0.5,9999.5);
 379
 380     hname="hNtrackletsOut";
 381     TH1F* hNtrackletsOut=new TH1F(hname.Data(),"Number of tracklets out of Centrality range;ntracklets;Entries",5000,-0.5,4999.5);
 382
 383     hname="hMultOut";
 384     TH1F* hMultOut=new TH1F(hname.Data(),"Multiplicity out of Centrality range;multiplicity;Entries",10000,-0.5,9999.5);
 385
 386     fOutputCheckCentrality->Add(hNtrackletsIn);
 387     fOutputCheckCentrality->Add(hNtrackletsOut);
 388     fOutputCheckCentrality->Add(hMultIn);
 389     fOutputCheckCentrality->Add(hMultOut);
 390
 391     PostData(6,fOutputCheckCentrality);
 392
 393   } else{
 394     hname="hNtracklets";
 395     TH1F* hNtracklets=new TH1F(hname.Data(),"Number of tracklets;ntracklets;Entries",5000,-0.5,4999.5);
 396
 397     hname="hMult";
 398     TH1F* hMult=new TH1F(hname.Data(),"Multiplicity;multiplicity;Entries",10000,-0.5,9999.5);
 399     fOutputTrack->Add(hNtracklets);
 400     fOutputTrack->Add(hMult);
 401   }
 402
 403   // Post the data
 404   PostData(1,fNEntries);
 405   PostData(2,fOutputPID);
 406   PostData(3,fOutputTrack);
 407   PostData(4,fCuts);
 408   PostData(5,fOutputCounters);
 409 }
 410
 411 //___________________________________________________________________________
 412 void AliAnalysisTaskSEHFQA::UserExec(Option_t */*option*/)
 413 {
 414   // Execute analysis for current event
 415
 416   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 417   if(fDebug>2) printf("Analysing decay %d\n",fDecayChannel);
 418   // Post the data already here
 419   PostData(1,fNEntries);
 420   PostData(2,fOutputPID);
 421   PostData(3,fOutputTrack);
 422   PostData(4,fCuts);
 423   PostData(5,fOutputCounters);
 424   if(fCuts->GetUseCentrality()) PostData(6,fOutputCheckCentrality);
 425
 426   TClonesArray *arrayProng =0;
 427   Int_t pdg=0;
 428   Int_t *pdgdaughters=0x0;
 429
 430   if(!aod && AODEvent() && IsStandardAOD()) {
 431     // In case there is an AOD handler writing a standard AOD, use the AOD
 432     // event in memory rather than the input (ESD) event.
 433     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 434     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 435     // have to taken from the AOD event hold by the AliAODExtension
 436     AliAODHandler* aodHandler = (AliAODHandler*)
 437       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 438     if(aodHandler->GetExtensions()) {
 439
 440       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 441       AliAODEvent *aodFromExt = ext->GetAOD();
 442
 443
 444
 445       switch(fDecayChannel){
 446       case 0:
 447         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 448         pdg=411;
 449         if(fReadMC){
 450           pdgdaughters =new Int_t[3];
 451           pdgdaughters[0]=211;//pi
 452           pdgdaughters[1]=321;//K
 453           pdgdaughters[2]=211;//pi
 454         }
 455         break;
 456       case 1:
 457         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 458         pdg=421;
 459         if(fReadMC){
 460           pdgdaughters =new Int_t[2];
 461           pdgdaughters[0]=211;//pi
 462           pdgdaughters[1]=321;//K
 463         }
 464         break;
 465       case 2:
 466         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 467         pdg=413;
 468         if(fReadMC){
 469           pdgdaughters =new Int_t[3];
 470           pdgdaughters[1]=211;//pi
 471           pdgdaughters[0]=321;//K
 472           pdgdaughters[2]=211;//pi (soft?)
 473         }
 474         break;
 475       case 3:
 476         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 477         pdg=431;
 478         if(fReadMC){
 479           pdgdaughters =new Int_t[3];
 480           pdgdaughters[0]=321;//K
 481           pdgdaughters[1]=321;//K
 482           pdgdaughters[2]=211;//pi
 483         }
 484         break;
 485       case 4:
 486         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 487         pdg=421;
 488         if(fReadMC){
 489           pdgdaughters =new Int_t[4];
 490           pdgdaughters[0]=321;
 491           pdgdaughters[1]=211;
 492           pdgdaughters[2]=211;
 493           pdgdaughters[3]=211;
 494         }
 495         break;
 496       case 5:
 497         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 498         pdg=4122;
 499         if(fReadMC){
 500           pdgdaughters =new Int_t[3];
 501           pdgdaughters[0]=2212;//p
 502           pdgdaughters[1]=321;//K
 503           pdgdaughters[2]=211;//pi
 504         }
 505         break;
 506       }
 507     }
 508   } else if(aod) {
 509     switch(fDecayChannel){
 510     case 0:
 511       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 512       pdg=411;
 513       if(fReadMC){
 514         pdgdaughters =new Int_t[3];
 515         pdgdaughters[0]=211;//pi
 516         pdgdaughters[1]=321;//K
 517         pdgdaughters[2]=211;//pi
 518       }
 519       break;
 520     case 1:
 521       arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 522       pdg=421;
 523       if(fReadMC){
 524         pdgdaughters =new Int_t[2];
 525         pdgdaughters[0]=211;//pi
 526         pdgdaughters[1]=321;//K
 527       }
 528       break;
 529     case 2:
 530       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
 531       pdg=413;
 532       if(fReadMC){
 533         pdgdaughters =new Int_t[3];
 534         pdgdaughters[1]=211;//pi
 535         pdgdaughters[0]=321;//K
 536         pdgdaughters[2]=211;//pi (soft?)
 537       }
 538       break;
 539     case 3:
 540       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 541       pdg=431;
 542       if(fReadMC){
 543         pdgdaughters =new Int_t[3];
 544         pdgdaughters[0]=321;//K
 545         pdgdaughters[1]=321;//K
 546         pdgdaughters[2]=211;//pi
 547       }
 548       break;
 549     case 4:
 550       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm4Prong");
 551       pdg=421;
 552       if(fReadMC){
 553         pdgdaughters =new Int_t[4];
 554         pdgdaughters[0]=321;
 555         pdgdaughters[1]=211;
 556         pdgdaughters[2]=211;
 557         pdgdaughters[3]=211;
 558       }
 559       break;
 560     case 5:
 561       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 562       pdg=4122;
 563       if(fReadMC){
 564         pdgdaughters =new Int_t[3];
 565         pdgdaughters[0]=2212;//p
 566         pdgdaughters[1]=321;//K
 567         pdgdaughters[2]=211;//pi
 568       }
 569       break;
 570     }
 571   }
 572
 573   if(!aod) {delete [] pdgdaughters;return;}
 574
 575
 576   Bool_t isSimpleMode=fSimpleMode;
 577   if(!arrayProng) {
 578     AliInfo("Branch not found! The output will contain only trak related histograms\n");
 579     isSimpleMode=kTRUE;
 580     fNEntries->Fill(2);
 581   }
 582
 583   TClonesArray *mcArray = 0;
 584   AliAODMCHeader *mcHeader = 0;
 585
 586   //check if MC
 587   if(fReadMC) {
 588     // load MC particles
 589     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 590     if(!mcArray) {
 591       printf("AliAnalysisTaskSEHFQA::UserExec: MC particles branch not found!\n");
 592       delete [] pdgdaughters;
 593       return;
 594     }
 595
 596     // load MC header
 597     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 598     if(!mcHeader) {
 599       printf("AliAnalysisTaskSEHFQA::UserExec: MC header branch not found!\n");
 600       delete [] pdgdaughters;
 601       return;
 602     }
 603   }
 604   // fix for temporary bug in ESDfilter
 605   // the AODs with null vertex pointer didn't pass the PhysSel
 606   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) {delete [] pdgdaughters; return;}
 607
 608   // count event
 609   fNEntries->Fill(0);
 610
 611   //count events with good vertex
 612   // AOD primary vertex
 613   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 614   TString primTitle = vtx1->GetTitle();
 615   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) fNEntries->Fill(4);
 616
 617   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
 618   TString trigclass=aod->GetFiredTriggerClasses();
 619   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNEntries->Fill(5);
 620
 621   Bool_t evSelbyCentrality=kTRUE,evSelected=kTRUE;
 622   //select event
 623   if(!fCuts->IsEventSelected(aod)) {
 624     evSelected=kFALSE;
 625     if(fCuts->GetWhyRejection()==1) fNEntries->Fill(1); // rejected for pileup
 626     if(fCuts->GetWhyRejection()==2) evSelbyCentrality=kFALSE; //rejected by centrality
 627   }
 628   if(evSelected || (!evSelected && !evSelbyCentrality)){ //events selected or not selected because of vtx
 629     if(fCuts->GetUseCentrality()){
 630       Int_t runNumber = aod->GetRunNumber();
 631       Int_t stdCent = (Int_t)(fCuts->GetCentrality(aod)+0.5);
 632       Int_t secondCent = (Int_t)(fCuts->GetCentrality(aod,fEstimator)+0.5);
 633       Int_t mincent=stdCent-stdCent%10;
 634       ((AliCounterCollection*)fOutputCounters->FindObject("stdEstimator"))->Count(Form("centralityclass:%d_%d/Run:%d",mincent,mincent+10,runNumber));
 635       mincent=secondCent-secondCent%10;
 636       ((AliCounterCollection*)fOutputCounters->FindObject("secondEstimator"))->Count(Form("centralityclass:%d_%d/Run:%d",mincent,mincent+10,runNumber));
 637
 638       if(stdCent<fCuts->GetMinCentrality() || stdCent>fCuts->GetMaxCentrality()){
 639         ((TH1F*)fOutputCheckCentrality->FindObject("hNtrackletsOut"))->Fill(aod->GetTracklets()->GetNumberOfTracklets());
 640         ((TH1F*)fOutputCheckCentrality->FindObject("hMultOut"))->Fill(aod->GetHeader()->GetRefMultiplicity());
 641       }else{
 642         ((TH1F*)fOutputCheckCentrality->FindObject("hNtrackletsIn"))->Fill(aod->GetTracklets()->GetNumberOfTracklets());
 643         ((TH1F*)fOutputCheckCentrality->FindObject("hMultIn"))->Fill(aod->GetHeader()->GetRefMultiplicity());
 644       }
 645
 646       PostData(6,fOutputCheckCentrality);
 647
 648     } else{
 649       ((TH1F*)fOutputTrack->FindObject("hNtracklets"))->Fill(aod->GetTracklets()->GetNumberOfTracklets());
 650       ((TH1F*)fOutputTrack->FindObject("hMult"))->Fill(aod->GetHeader()->GetRefMultiplicity());
 651     }
 652   }
 653
 654   if(!evSelected) {
 655     delete [] pdgdaughters;
 656     return; //discard all events not selected (vtx and/or centrality)
 657   }
 658
 659   Int_t ntracks=0;
 660   Int_t isGoodTrack=0;
 661
 662   if(aod) ntracks=aod->GetNTracks();
 663
 664   //loop on tracks in the event
 665   for (Int_t k=0;k<ntracks;k++){
 666     AliAODTrack* track=aod->GetTrack(k);
 667     AliAODPidHF* pidHF=fCuts->GetPidHF();
 668     AliAODPid *pid = track->GetDetPid();
 669     if(!pid)  {if (fDebug>1)cout<<"No AliAODPid found"<<endl; continue;}
 670
 671     Double_t times[AliPID::kSPECIES];
 672     pid->GetIntegratedTimes(times);
 673
 674     Double_t tofRes[AliPID::kSPECIES];
 675     pid->GetTOFpidResolution(tofRes);
 676
 677     //check TOF
 678     if(pidHF && pidHF->CheckStatus(track,"TOF")){
 679       ((TH1F*)fOutputPID->FindObject("hTOFtime"))->Fill(times[AliPID::kProton]);
 680       ((TH2F*)fOutputPID->FindObject("hTOFtimeKaonHyptime"))->Fill(track->P(),pid->GetTOFsignal()-times[3]); //3 is kaon
 681       ((TH1F*)fOutputPID->FindObject("hTOFsig"))->Fill(pid->GetTOFsignal());
 682       if (pid->GetTOFsignal()< 0) ((TH1F*)fOutputPID->FindObject("hTOFsig"))->Fill(-1);
 683
 684       // test a "simple" 160 ps TOF sigma PID
 685       ((TH2F*)fOutputPID->FindObject("hTOFsigmaK160"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kKaon])/160);
 686       ((TH2F*)fOutputPID->FindObject("hTOFsigmaPion160"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kPion])/160);
 687       ((TH2F*)fOutputPID->FindObject("hTOFsigmaProton160"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kProton])/160);
 688
 689       // test TOF sigma PID
 690       if (tofRes[2] != 0.) {   // protection against 'old' AODs...
 691         ((TH2F*)fOutputPID->FindObject("hTOFsigmaKSigPid"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kKaon])/tofRes[3]);
 692         ((TH2F*)fOutputPID->FindObject("hTOFsigmaPionSigPid"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kPion])/tofRes[2]);
 693         ((TH2F*)fOutputPID->FindObject("hTOFsigmaProtonSigPid"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kProton])/tofRes[4]);
 694         for (Int_t iS=2; iS<5; iS++){ //we plot TOF Pid resolution for 3-sigma identified particles
 695           if ( (TMath::Abs(times[iS]-pid->GetTOFsignal())/tofRes[iS])<3.){
 696             switch (iS) {
 697             case AliPID::kPion:
 698               ((TH1F*)fOutputPID->FindObject("hTOFsigPid3sigPion"))->Fill(tofRes[iS]);
 699               break;
 700             case AliPID::kKaon:
 701               ((TH1F*)fOutputPID->FindObject("hTOFsigPid3sigKaon"))->Fill(tofRes[iS]);
 702               break;
 703             case AliPID::kProton:
 704               ((TH1F*)fOutputPID->FindObject("hTOFsigPid3sigProton"))->Fill(tofRes[iS]);
 705               break;
 706             default:
 707               break;
 708             }
 709           }
 710         }
 711       }
 712     }//if TOF status
 713
 714     if(pidHF && pidHF->CheckStatus(track,"TPC")){
 715       Double_t alephParameters[5];
 716       //this is recommended for LHC10d
 717       alephParameters[0] = 1.34490e+00/50.;
 718       alephParameters[1] =  2.69455e+01;
 719       alephParameters[2] =  TMath::Exp(-2.97552e+01);
 720       alephParameters[3] = 2.35339e+00;
 721       alephParameters[4] = 5.98079e+00;
 722       /*
 723       //this is recommended for LHC10bc
 724       alephParameters[0] = 0.0283086/0.97;
 725       alephParameters[1] = 2.63394e+01;
 726       alephParameters[2] = 5.04114e-11;
 727       alephParameters[3] = 2.12543e+00;
 728       alephParameters[4] = 4.88663e+00;
 729       */
 730       AliTPCPIDResponse* tpcres=new AliTPCPIDResponse();
 731       tpcres->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
 732       Double_t TPCp=pid->GetTPCmomentum();
 733       Double_t TPCsignal=pid->GetTPCsignal();
 734       ((TH1F*)fOutputPID->FindObject("hTPCsig"))->Fill(TPCsignal);
 735       ((TH1F*)fOutputPID->FindObject("hTPCsigvsp"))->Fill(TPCp,TPCsignal);
 736       //if (pidHF->IsKaonRaw(track, "TOF"))
 737       ((TH2F*)fOutputPID->FindObject("hTPCsigmaK"))->Fill(TPCp,tpcres->GetNumberOfSigmas(TPCp,TPCsignal,track->GetTPCNcls(),AliPID::kKaon));
 738       //if (pidHF->IsPionRaw(track, "TOF"))
 739       ((TH2F*)fOutputPID->FindObject("hTPCsigmaPion"))->Fill(TPCp,tpcres->GetNumberOfSigmas(TPCp,TPCsignal,track->GetTPCNcls(),AliPID::kPion));
 740       //if (pidHF->IsProtonRaw(track,"TOF"))
 741       ((TH2F*)fOutputPID->FindObject("hTPCsigmaProton"))->Fill(TPCp,tpcres->GetNumberOfSigmas(TPCp,TPCsignal,track->GetTPCNcls(),AliPID::kProton));
 742       delete tpcres;
 743     }//if TPC status
 744
 745     //check clusters of the tracks
 746     Int_t nclsTot=0,nclsSPD=0;
 747
 748     for(Int_t l=0;l<6;l++) {
 749       if(TESTBIT(track->GetITSClusterMap(),l)) {
 750         nclsTot++; if(l<2) nclsSPD++;
 751       }
 752     }
 753     ((TH1F*)fOutputTrack->FindObject("hnClsITS"))->Fill(nclsTot);
 754     ((TH1F*)fOutputTrack->FindObject("hnClsSPD"))->Fill(nclsSPD);
 755
 756     if(!(track->GetStatus()&AliESDtrack::kTPCin) && track->GetStatus()&AliESDtrack::kITSrefit && !(track->GetStatus()&AliESDtrack::kITSpureSA)){//tracks retrieved in the ITS and not reconstructed in the TPC
 757       ((TH1F*)fOutputTrack->FindObject("hnClsITS-SA"))->Fill(nclsTot);
 758     }
 759
 760     if(isSimpleMode){
 761
 762       if (track->Pt()>0.3 &&
 763           track->GetStatus()&AliESDtrack::kTPCrefit &&
 764           track->GetStatus()&AliESDtrack::kITSrefit &&
 765           /*nclsTot>3 &&*/
 766           nclsSPD>0) {//fill hist good tracks
 767
 768         ((TH1F*)fOutputTrack->FindObject("hptGoodTr"))->Fill(track->Pt());
 769
 770         isGoodTrack++;
 771
 772         ((TH1F*)fOutputTrack->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
 773       }
 774     }//simple mode: no IsSelected on tracks: use "manual" cuts
 775
 776   } //end loop on tracks
 777
 778   if(!isSimpleMode){
 779     // loop over candidates
 780     Int_t nCand = arrayProng->GetEntriesFast();
 781     Int_t ndaugh=3;
 782     if(fDecayChannel==AliAnalysisTaskSEHFQA::kD0toKpi) ndaugh=2;
 783     if(fDecayChannel==AliAnalysisTaskSEHFQA::kD0toKpipipi) ndaugh=4;
 784
 785     for (Int_t iCand = 0; iCand < nCand; iCand++) {
 786       AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayProng->UncheckedAt(iCand);
 787
 788       if(fReadMC){
 789         Int_t labD = d->MatchToMC(pdg,mcArray,ndaugh,pdgdaughters);
 790         if(labD>=0){
 791           AliAODMCParticle *partD = (AliAODMCParticle*)mcArray->At(labD);
 792           Int_t label=partD->GetMother();
 793           AliAODMCParticle *mot = (AliAODMCParticle*)mcArray->At(label);
 794           while(label>=0){//get first mother
 795             mot = (AliAODMCParticle*)mcArray->At(label);
 796             label=mot->GetMother();
 797           }
 798           Int_t pdgMotCode = mot->GetPdgCode();
 799
 800           if(TMath::Abs(pdgMotCode)==4) fNEntries->Fill(6); //from primary charm
 801           if(TMath::Abs(pdgMotCode)==5) fNEntries->Fill(7); //from beauty
 802
 803         }
 804       }//end MC
 805       else fNEntries->Fill(6); //count the candidates (data)
 806
 807       for(Int_t id=0;id<ndaugh;id++){
 808
 809         //other histograms to be filled when the cut object is given
 810         AliAODTrack* track=(AliAODTrack*)d->GetDaughter(id);
 811         if(fReadMC){
 812           Int_t label=track->GetLabel();
 813           if (label<0)fNEntries->Fill(8);
 814         }
 815         //track quality
 816
 817         if (fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(pdg)) && fCuts->IsSelected(d,AliRDHFCuts::kTracks,aod)) {
 818
 819           ((TH1F*)fOutputTrack->FindObject("hptGoodTr"))->Fill(track->Pt());
 820           isGoodTrack++;
 821
 822           ((TH1F*)fOutputTrack->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
 823
 824           ((TH1F*)fOutputTrack->FindObject("hd0"))->Fill(d->Getd0Prong(id));
 825
 826           if (fCuts->IsSelected(d,AliRDHFCuts::kAll,aod)){
 827
 828             AliAODPid *pid = track->GetDetPid();
 829             AliAODPidHF* pidHF=fCuts->GetPidHF();
 830             Double_t times[5];
 831             pid->GetIntegratedTimes(times);
 832             if(pidHF && pidHF->CheckStatus(track,"TOF")) ((TH2F*)fOutputPID->FindObject("hTOFtimeKaonHyptimeAC"))->Fill(track->P(),pid->GetTOFsignal()-times[AliPID::kKaon]);
 833             if(pidHF && pidHF->CheckStatus(track,"TPC")) ((TH2F*)fOutputPID->FindObject("hTPCsigvspAC"))->Fill(pid->GetTPCmomentum(),pid->GetTPCsignal());
 834
 835             fNEntries->Fill(3);
 836           } //end analysis cuts
 837         } //end acceptance and track cuts
 838       } //end loop on tracks in the candidate
 839     } //end loop on candidates
 840   }
 841
 842   delete [] pdgdaughters;
 843   PostData(1,fNEntries);
 844   PostData(2,fOutputPID);
 845   PostData(3,fOutputTrack);
 846   PostData(4,fCuts);
 847   PostData(5,fOutputCounters);
 848   //Post data 6 done in case of centrality on
 849
 850 }
 851
 852 //____________________________________________________________________________
 853 void AliAnalysisTaskSEHFQA::Terminate(Option_t */*option*/){
 854   //terminate analysis
 855
 856   fNEntries = dynamic_cast<TH1F*>(GetOutputData(1));
 857   if(!fNEntries){
 858     printf("ERROR: %s not available\n",GetOutputSlot(1)->GetContainer()->GetName());
 859     return;
 860   }
 861
 862   fOutputPID = dynamic_cast<TList*> (GetOutputData(2));
 863   if (!fOutputPID) {
 864     printf("ERROR: %s not available\n",GetOutputSlot(2)->GetContainer()->GetName());
 865     return;
 866   }
 867
 868   fOutputTrack = dynamic_cast<TList*> (GetOutputData(3));
 869   if (!fOutputTrack) {
 870     printf("ERROR: %s not available\n",GetOutputSlot(3)->GetContainer()->GetName());
 871     return;
 872   }
 873
 874 }
 875