PWGHF/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 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // AliAnalysisTaskSE for HF quality assurance
  21 //
  22 // Author: Chiara Bianchin, chiara.bianchin@pd.infn.it
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <Riostream.h>
  26 #include <TClonesArray.h>
  27 #include <TCanvas.h>
  28 #include <TList.h>
  29 #include <TH1F.h>
  30 #include <TH2F.h>
  31 #include <TH3F.h>
  32 #include <TProfile2D.h>
  33 #include <TDatabasePDG.h>
  34
  35 #include <AliAnalysisDataSlot.h>
  36 #include <AliAnalysisDataContainer.h>
  37 #include "AliAnalysisManager.h"
  38 #include "AliESDtrack.h"
  39 #include "AliESDVertex.h"
  40 #include "AliVertexerTracks.h"
  41 #include "AliPID.h"
  42 #include "AliTPCPIDResponse.h"
  43 #include "AliAODHandler.h"
  44 #include "AliAODEvent.h"
  45 #include "AliAODVertex.h"
  46 #include "AliAODTrack.h"
  47 #include "AliAODMCParticle.h"
  48 #include "AliAODMCHeader.h"
  49 #include "AliAODRecoDecayHF2Prong.h"
  50 #include "AliAODRecoCascadeHF.h"
  51 #include "AliAnalysisVertexingHF.h"
  52 #include "AliAnalysisTaskSE.h"
  53 #include "AliCounterCollection.h"
  54 #include "AliRDHFCuts.h"
  55 #include "AliRDHFCutsDplustoKpipi.h"
  56 #include "AliRDHFCutsD0toKpipipi.h"
  57 #include "AliRDHFCutsDstoKKpi.h"
  58 #include "AliRDHFCutsDStartoKpipi.h"
  59 #include "AliRDHFCutsD0toKpi.h"
  60 #include "AliRDHFCutsLctopKpi.h"
  61 #include "AliRDHFCutsLctoV0.h"
  62 #include "AliInputEventHandler.h"
  63
  64 #include "AliFlowEvent.h"
  65 #include "AliFlowTrackCuts.h"
  66 #include "AliFlowTrackSimple.h"
  67 #include "AliFlowVector.h"
  68
  69 #include "AliAnalysisTaskSEHFQA.h"
  70
  71 using std::cout;
  72 using std::endl;
  73
  74 ClassImp(AliAnalysisTaskSEHFQA)
  75
  76 //____________________________________________________________________________
  77
  78 AliAnalysisTaskSEHFQA::AliAnalysisTaskSEHFQA():AliAnalysisTaskSE(),
  79   fNEntries(0x0),
  80   fOutputPID(0x0),
  81   fOutputTrack(0x0),
  82   fOutputCounters(0x0),
  83   fOutputCheckCentrality(0x0),
  84   fOutputEvSelection(0x0),
  85   fOutputFlowObs(0x0),
  86   fDecayChannel(AliAnalysisTaskSEHFQA::kD0toKpi),
  87   fCuts(0x0),
  88   fFlowEvent(0x0),
  89   fRFPcuts(0x0),
  90   fEstimator(AliRDHFCuts::kCentTRK),
  91   fReadMC(kFALSE),
  92   fSimpleMode(kFALSE),
  93   fOnOff()
  94 {
  95   //default constructor
  96   fOnOff[0]=kTRUE;
  97   fOnOff[1]=kTRUE;
  98   fOnOff[2]=kTRUE;
  99   fOnOff[3]=kTRUE;
 100   fOnOff[4]=kTRUE;
 101 }
 102
 103 //____________________________________________________________________________
 104 AliAnalysisTaskSEHFQA::AliAnalysisTaskSEHFQA(const char *name, AliAnalysisTaskSEHFQA::DecChannel ch,AliRDHFCuts* cuts):
 105   AliAnalysisTaskSE(name),
 106   fNEntries(0x0),
 107   fOutputPID(0x0),
 108   fOutputTrack(0x0),
 109   fOutputCounters(0x0),
 110   fOutputCheckCentrality(0x0),
 111   fOutputEvSelection(0x0),
 112   fOutputFlowObs(0x0),
 113   fDecayChannel(ch),
 114   fCuts(0x0),
 115   fFlowEvent(0x0),
 116   fRFPcuts(0x0),
 117   fEstimator(AliRDHFCuts::kCentTRK),
 118   fReadMC(kFALSE),
 119   fSimpleMode(kFALSE),
 120   fOnOff()
 121 {
 122   //constructor
 123
 124   //SetCutObject(cuts);
 125   fCuts=cuts;
 126
 127   fOnOff[0]=kTRUE;
 128   fOnOff[1]=kTRUE;
 129   fOnOff[2]=kTRUE;
 130   fOnOff[3]=kTRUE;
 131   fOnOff[4]=kTRUE;
 132
 133   // Output slot #1 writes into a TH1F container (number of events)
 134   DefineOutput(1,TH1F::Class());  //My private output
 135   // Output slot #2 writes into a TList container (PID)
 136   if (fOnOff[1]) DefineOutput(2,TList::Class());  //My private output
 137   // Output slot #3 writes into a TList container (Tracks)
 138   if (fOnOff[0]) DefineOutput(3,TList::Class());  //My private output
 139   // Output slot #4 writes into a AliRDHFCuts container (cuts)
 140   switch(fDecayChannel){
 141   case 0:
 142     DefineOutput(4,AliRDHFCutsDplustoKpipi::Class());  //My private output
 143     break;
 144   case 1:
 145     DefineOutput(4,AliRDHFCutsD0toKpi::Class());  //My private output
 146     break;
 147   case 2:
 148     DefineOutput(4,AliRDHFCutsDStartoKpipi::Class());  //My private output
 149     break;
 150   case 3:
 151     DefineOutput(4,AliRDHFCutsDstoKKpi::Class());  //My private output
 152     break;
 153   case 4:
 154     DefineOutput(4,AliRDHFCutsD0toKpipipi::Class());  //My private output
 155     break;
 156   case 5:
 157     DefineOutput(4,AliRDHFCutsLctopKpi::Class());  //My private output
 158     break;
 159   case kLambdactoV0:
 160     DefineOutput(4,AliRDHFCutsLctoV0::Class());  //My private output
 161     break;
 162   }
 163   if (fOnOff[2]) {
 164     // Output slot #5 writes into a TList container (AliCounterCollection)
 165     DefineOutput(5,TList::Class());  //My private output
 166     // Output slot #6 writes into a TList container (TH1F)
 167     DefineOutput(6,TList::Class());  //My private output
 168   }
 169
 170   if(fOnOff[3]) DefineOutput(7,TList::Class());  //My private output
 171   if(fOnOff[4]) DefineOutput(8,TList::Class());  //My private output
 172
 173 }
 174
 175 //___________________________________________________________________________
 176 AliAnalysisTaskSEHFQA::~AliAnalysisTaskSEHFQA()
 177 {
 178   //destructor
 179   delete fNEntries;
 180
 181   delete fOutputPID;
 182
 183   delete fOutputTrack;
 184
 185   delete fOutputCounters;
 186
 187   delete fOutputCheckCentrality;
 188
 189   delete fOutputEvSelection;
 190
 191   if(fOnOff[4]) {
 192     delete fOutputFlowObs;
 193     delete fFlowEvent;
 194   }
 195 }
 196
 197 //___________________________________________________________________________
 198 void AliAnalysisTaskSEHFQA::Init(){
 199
 200   //initialization
 201   if(fDebug > 1) printf("AnalysisTaskSEHFQA::Init() \n");
 202
 203   switch(fDecayChannel){
 204   case 0:
 205     {
 206       AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
 207       // Post the data
 208       PostData(4,copycut);
 209     }
 210     break;
 211   case 1:
 212     {
 213       AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 214       // Post the data
 215       PostData(4,copycut);
 216     }
 217     break;
 218   case 2:
 219     {
 220       AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 221       // Post the data
 222       PostData(4,copycut);
 223     }
 224     break;
 225   case 3:
 226     {
 227       AliRDHFCutsDstoKKpi* copycut=new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
 228       // Post the data
 229       PostData(4,copycut);
 230     }
 231     break;
 232   case 4:
 233     {
 234       AliRDHFCutsD0toKpipipi* copycut=new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
 235       // Post the data
 236       PostData(4,copycut);
 237     }
 238     break;
 239   case 5:
 240     {
 241       AliRDHFCutsLctopKpi* copycut=new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
 242       // Post the data
 243       PostData(4,copycut);
 244     }
 245     break;
 246   case kLambdactoV0:
 247     {
 248       AliRDHFCutsLctoV0* copycut=new AliRDHFCutsLctoV0(*(static_cast<AliRDHFCutsLctoV0*>(fCuts)));
 249       // Post the data
 250       PostData(4,copycut);
 251     }
 252     break;
 253
 254   default:
 255     return;
 256   }
 257
 258
 259
 260 }
 261
 262 //___________________________________________________________________________
 263 void AliAnalysisTaskSEHFQA::UserCreateOutputObjects()
 264 {
 265
 266   //create the output container
 267   if(fDebug > 1) printf("AnalysisTaskSEHFQA::UserCreateOutputObjects() \n");
 268
 269   //count events
 270
 271   fNEntries=new TH1F(GetOutputSlot(1)->GetContainer()->GetName(), "Counts the number of events", 10,-0.5,9.5);
 272   fNEntries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 273   fNEntries->GetXaxis()->SetBinLabel(2,"Pile-up Rej");
 274   fNEntries->GetXaxis()->SetBinLabel(3,"No VertexingHF");
 275   fNEntries->GetXaxis()->SetBinLabel(4,"nCandidates(AnCuts)");
 276   fNEntries->GetXaxis()->SetBinLabel(5,"EventsWithGoodVtx");
 277   //fNEntries->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
 278   fNEntries->GetXaxis()->SetBinLabel(6,"N. of CSH1");
 279   if(fReadMC){
 280     fNEntries->GetXaxis()->SetBinLabel(7,"MC Cand from c");
 281     fNEntries->GetXaxis()->SetBinLabel(8,"MC Cand from b");
 282     fNEntries->GetXaxis()->SetBinLabel(9,"N fake Trks");
 283     fNEntries->GetXaxis()->SetBinLabel(10,"N true Trks");
 284   } else{
 285     fNEntries->GetXaxis()->SetBinLabel(7,"N candidates");
 286   }
 287
 288   fNEntries->GetXaxis()->SetNdivisions(1,kFALSE);
 289
 290   //PID
 291   if(fOnOff[1]){
 292     fOutputPID=new TList();
 293     fOutputPID->SetOwner();
 294     fOutputPID->SetName(GetOutputSlot(2)->GetContainer()->GetName());
 295
 296     //TOF pid
 297     TH1F* hTOFflags=new TH1F("hTOFflags","TOF flags",6,-0.5,5.5);
 298     hTOFflags->SetMinimum(0.);
 299     hTOFflags->GetXaxis()->SetBinLabel(1,"All Tracks");
 300     hTOFflags->GetXaxis()->SetBinLabel(2,"kTPCout");
 301     hTOFflags->GetXaxis()->SetBinLabel(3,"kTOFout");
 302     hTOFflags->GetXaxis()->SetBinLabel(4,"kTIME");
 303     hTOFflags->GetXaxis()->SetBinLabel(5,"kTOFpid");
 304     hTOFflags->GetXaxis()->SetBinLabel(6,"kTOFmismatch");
 305
 306     TString hname="hTOFsig";
 307     TH1F* hTOFsig=new TH1F(hname.Data(),"Distribution of TOF signal;TOF time [ps];Entries", 100, -2.e3,40.e3);
 308
 309     hname="hTOFtime";
 310     TH1F* hTOFtime=new TH1F(hname.Data(),"Distribution of TOF time Kaon;TOF time(Kaon) [ps];Entries", 1000, 0.,50000.);
 311
 312     hname="hTOFtimeKaonHyptime";
 313     TH2F* hTOFtimeKaonHyptime=new TH2F(hname.Data(),"TOFtime - timeHypothesisForKaon;p[GeV/c];TOFtime - timeHypothesisForKaon [ps]",500,0.,10.,1000,-20000.,20000.);
 314
 315     hname="hTOFtimeKaonHyptimeAC";
 316     TH2F* hTOFtimeKaonHyptimeAC=new TH2F(hname.Data(),"TOFtime - timeHypothesisForKaon;p[GeV/c];TOFtime - timeHypothesisForKaon [ps]",500,0.,10.,1000,-20000.,20000.);
 317
 318     hname="hTOFsigmaKSigPid";
 319     TH2F* hTOFsigmaKSigPid=new TH2F(hname.Data(),"(TOFsignal-timeK)/tofSigPid;p[GeV/c];(TOFsignal-timeK)/tofSigPid",500,0.,10.,400,-20,20);
 320
 321     hname="hTOFsigmaPionSigPid";
 322     TH2F* hTOFsigmaPionSigPid=new TH2F(hname.Data(),"(TOFsignal-time#pi)/tofSigPid;p[GeV/c];(TOFsignal-time#pi)/tofSigPid",500,0.,10.,400,-20,20);
 323
 324     hname="hTOFsigmaProtonSigPid";
 325     TH2F* hTOFsigmaProtonSigPid=new TH2F(hname.Data(),"(TOFsignal-timep)/tofSigPid;p[GeV/c];(TOFsignal-time p)/tofSigPid",500,0.,10.,400,-20,20);
 326
 327     hname="hTOFsigPid3sigPion";
 328     TH1F* hTOFsigPid3sigPion=new TH1F(hname.Data(),"TOF PID resolution (#pi) [ps]",500,0.,1000.);
 329
 330     hname="hTOFsigPid3sigKaon";
 331     TH1F* hTOFsigPid3sigKaon=new TH1F(hname.Data(),"TOF PID resolution (K) [ps]",500,0.,1000.);
 332
 333     hname="hTOFsigPid3sigProton";
 334     TH1F* hTOFsigPid3sigProton=new TH1F(hname.Data(),"TOF PID resolution (p) [ps]",500,0.,1000.);
 335
 336
 337     //TPC pid
 338     hname="hTPCsig";
 339     TH1F* hTPCsig=new TH1F(hname.Data(),"Distribution of TPC signal;TPC sig;Entries", 100, 35.,100.);
 340
 341     hname="hTPCsigvsp";
 342     TH2F* hTPCsigvsp=new TH2F(hname.Data(),"TPCsig vs p;TPC p[GeV/c];TPCsig",500,0.,10.,1000,35.,100.);
 343
 344     hname="hTPCsigvspAC";
 345     TH2F* hTPCsigvspAC=new TH2F(hname.Data(),"TPCsig vs p;TPCp[GeV/c];TPCsig",500,0.,10.,1000,35.,100.);
 346
 347     hname="hTPCsigmaK";
 348     TH2F* hTPCsigmaK=new TH2F(hname.Data(),"TPC Sigma for K as a function of momentum;p[GeV/c];Sigma Kaon",500,0.,10.,400,-20,20);
 349
 350     hname="hTPCsigmaPion";
 351     TH2F* hTPCsigmaPion=new TH2F(hname.Data(),"TPC Sigma for #pi as a function of momentum;p[GeV/c];Sigma #pi",500,0.,10.,400,-20,20);
 352
 353     hname="hTPCsigmaProton";
 354     TH2F* hTPCsigmaProton=new TH2F(hname.Data(),"TPC Sigma for proton as a function of momentum;p[GeV/c];Sigma Proton",500,0.,10.,400,-20,20);
 355
 356     fOutputPID->Add(hTOFflags);
 357     fOutputPID->Add(hTOFsig);
 358     fOutputPID->Add(hTPCsig);
 359     fOutputPID->Add(hTOFtime);
 360     fOutputPID->Add(hTOFtimeKaonHyptime);
 361     fOutputPID->Add(hTOFtimeKaonHyptimeAC);
 362     fOutputPID->Add(hTOFsigmaKSigPid);
 363     fOutputPID->Add(hTOFsigmaPionSigPid);
 364     fOutputPID->Add(hTOFsigmaProtonSigPid);
 365     fOutputPID->Add(hTOFsigPid3sigPion);
 366     fOutputPID->Add(hTOFsigPid3sigKaon);
 367     fOutputPID->Add(hTOFsigPid3sigProton);
 368     fOutputPID->Add(hTPCsigvsp);
 369     fOutputPID->Add(hTPCsigvspAC);
 370     fOutputPID->Add(hTPCsigmaK);
 371     fOutputPID->Add(hTPCsigmaPion);
 372     fOutputPID->Add(hTPCsigmaProton);
 373
 374     if(fReadMC){
 375       //TOF
 376       hname="hTOFsigmaMCKSigPid";
 377       TH2F* hTOFsigmaMCKSigPid=new TH2F(hname.Data(),"(TOFsignal-timeK)/tofSigPid;p[GeV/c];(TOFsignal-timeK)/tofSigPid",500,0.,10.,400,-20,20);
 378
 379       hname="hTOFsigmaMCPionSigPid";
 380       TH2F* hTOFsigmaMCPionSigPid=new TH2F(hname.Data(),"(TOFsignal-time#pi)/tofSigPid;p[GeV/c];(TOFsignal-time#pi)/tofSigPid",500,0.,10.,400,-20,20);
 381
 382       hname="hTOFsigmaMCProtonSigPid";
 383       TH2F* hTOFsigmaMCProtonSigPid=new TH2F(hname.Data(),"(TOFsignal-timep)/tofSigPid;p[GeV/c];(TOFsignal-time p)/tofSigPid",500,0.,10.,400,-20,20);
 384
 385       //TPC
 386       hname="hTPCsigmaMCK";
 387       TH2F* hTPCsigmaMCK=new TH2F(hname.Data(),"TPC Sigma for K as a function of momentum;p[GeV/c];Sigma Kaon",500,0.,10.,400,-20,20);
 388
 389       hname="hTPCsigmaMCPion";
 390       TH2F* hTPCsigmaMCPion=new TH2F(hname.Data(),"TPC Sigma for #pi as a function of momentum;p[GeV/c];Sigma #pi",500,0.,10.,400,-20,20);
 391
 392       hname="hTPCsigmaMCProton";
 393       TH2F* hTPCsigmaMCProton=new TH2F(hname.Data(),"TPC Sigma for proton as a function of momentum;p[GeV/c];Sigma Proton",500,0.,10.,400,-20,20);
 394
 395       fOutputPID->Add(hTOFsigmaMCKSigPid);
 396       fOutputPID->Add(hTOFsigmaMCPionSigPid);
 397       fOutputPID->Add(hTOFsigmaMCProtonSigPid);
 398       fOutputPID->Add(hTPCsigmaMCK);
 399       fOutputPID->Add(hTPCsigmaMCPion);
 400       fOutputPID->Add(hTPCsigmaMCProton);
 401
 402     }
 403   }
 404
 405   //quality of the tracks
 406   if(fOnOff[0]){
 407     fOutputTrack=new TList();
 408     fOutputTrack->SetOwner();
 409     fOutputTrack->SetName(GetOutputSlot(3)->GetContainer()->GetName());
 410
 411     TString hname="hnClsITS";
 412     TH1F* hnClsITS=new TH1F(hname.Data(),"Distribution of number of ITS clusters;nITScls;Entries",7,-0.5,6.5);
 413
 414     hname="hnClsITSselTr";
 415     TH1F* hnClsITSselTr=new TH1F(hname.Data(),"Distribution of number of ITS clusters selected tracks;nITScls;Entries",7,-0.5,6.5);
 416
 417     hname="hnClsITS-SA";
 418     TH1F* hnClsITSSA=new TH1F(hname.Data(),"Distribution of number of ITS clusters(ITS-SA);nITScls;Entries",7,-0.5,6.5);
 419
 420
 421     hname="hnLayerITS";
 422     TH1F* hnLayerITS=new TH1F(hname.Data(),"Number of tracks with point in layer;ITS layer;",7,-1.5,5.5);
 423     hnLayerITS->GetXaxis()->SetBinLabel(1,"n tracks");
 424     hnLayerITS->GetXaxis()->SetBinLabel(2,"SPDin");
 425     hnLayerITS->GetXaxis()->SetBinLabel(3,"SPDout");
 426     hnLayerITS->GetXaxis()->SetBinLabel(4,"SDDin");
 427     hnLayerITS->GetXaxis()->SetBinLabel(5,"SDDout");
 428     hnLayerITS->GetXaxis()->SetBinLabel(6,"SSDin");
 429     hnLayerITS->GetXaxis()->SetBinLabel(7,"SSDout");
 430
 431     hname="hnLayerITSsa";
 432     TH1F* hnLayerITSsa=new TH1F(hname.Data(),"Number of tracks with point in layer;ITS layer;",7,-1.5,5.5);
 433     hnLayerITSsa->GetXaxis()->SetBinLabel(1,"n tracks");
 434     hnLayerITSsa->GetXaxis()->SetBinLabel(2,"SPDin");
 435     hnLayerITSsa->GetXaxis()->SetBinLabel(3,"SPDout");
 436     hnLayerITSsa->GetXaxis()->SetBinLabel(4,"SDDin");
 437     hnLayerITSsa->GetXaxis()->SetBinLabel(5,"SDDout");
 438     hnLayerITSsa->GetXaxis()->SetBinLabel(6,"SSDin");
 439     hnLayerITSsa->GetXaxis()->SetBinLabel(7,"SSDout");
 440
 441     hname="hnClsSPD";
 442     TH1F* hnClsSPD=new TH1F(hname.Data(),"Distribution of number of SPD clusters;nSPDcls;Entries",3,-0.5,2.5);
 443
 444     hname="hptGoodTr";
 445     TH1F* hptGoodTr=new TH1F(hname.Data(),"Pt distribution of 'good' tracks;p_{t}[GeV];Entries/0.05 GeV/c",400,0.,20.);
 446     hptGoodTr->SetTitleOffset(1.3,"Y");
 447
 448     if(!fSimpleMode){
 449       hname="hptGoodTrFromDaugh";
 450       TH1F* hptGoodTrFromDaugh=new TH1F(hname.Data(),"Pt distribution of 'good' candidate's daughters;p_{t}[GeV];Entries/0.05 GeV/c",400,0.,20.);
 451       hptGoodTrFromDaugh->SetTitleOffset(1.3,"Y");
 452       fOutputTrack->Add(hptGoodTrFromDaugh);
 453     }
 454
 455     hname="hdistrGoodTr";
 456     TH1F* hdistrGoodTr=new TH1F(hname.Data(),"Distribution of number of 'good' candidate's daughters per event;no.good-tracks/ev;Entries",4000,-0.5,3999.5);
 457     hdistrGoodTr->SetTitleOffset(1.3,"Y");
 458
 459     hname="hdistrSelTr";
 460     TH1F* hdistrSelTr=new TH1F(hname.Data(),"Distribution of number of Selected tracks per event;no.good-tracks/ev;Entries",4000,-0.5,3999.5);
 461     hdistrSelTr->SetTitleOffset(1.3,"Y");
 462
 463     hname="hd0";
 464     TH1F* hd0=new TH1F(hname.Data(),"Impact parameter (rphi) distribution of 'good' tracks;d_{0rphi}[cm];Entries/10^{3} cm",200,-0.1,0.1);
 465
 466     hname="hd0z";
 467     TH1F* hd0z=new TH1F(hname.Data(),"Impact parameter (z) distribution of 'good' tracks;d_{0z}[cm];Entries/10^{3} cm",200,-0.1,0.1);
 468
 469     fOutputTrack->Add(hnClsITS);
 470     fOutputTrack->Add(hnClsITSselTr);
 471     fOutputTrack->Add(hnClsITSSA);
 472     fOutputTrack->Add(hnLayerITS);
 473     fOutputTrack->Add(hnLayerITSsa);
 474     fOutputTrack->Add(hnClsSPD);
 475     fOutputTrack->Add(hptGoodTr);
 476     fOutputTrack->Add(hdistrGoodTr);
 477     fOutputTrack->Add(hdistrSelTr);
 478     fOutputTrack->Add(hd0);
 479     fOutputTrack->Add(hd0z);
 480
 481     if(fReadMC){
 482       hname="hdistrFakeTr";
 483       TH1F* hdistrFakeTr=new TH1F(hname.Data(),"Distribution of number of fake tracks per event;no.fake-tracks/ev;Entries",4000,-0.5,3999.5);
 484       hdistrFakeTr->SetTitleOffset(1.3,"Y");
 485
 486       hname="hd0f";
 487       TH1F* hd0f=new TH1F(hname.Data(),"Impact parameter distribution of fake tracks;d_{0}[cm];Entries/10^{3} cm",200,-0.1,0.1);
 488
 489       hname="hptFakeTr";
 490       TH1F* hptFakeTr=new TH1F(hname.Data(),"Pt distribution of fake tracks;p_{t}[GeV];Entries/0.05 GeV/c",400,0.,20.);
 491       hptFakeTr->SetTitleOffset(1.3,"Y");
 492       if(!fSimpleMode){
 493         hname="hptFakeTrFromDaugh";
 494         TH1F* hptFakeTrFromDaugh=new TH1F(hname.Data(),"Pt distribution of fake tracks from daughters;p_{t}[GeV];Entries/0.05 GeV/c",400,0.,20.);
 495         hptFakeTrFromDaugh->SetTitleOffset(1.3,"Y");
 496         fOutputTrack->Add(hptFakeTrFromDaugh);
 497       }
 498
 499       fOutputTrack->Add(hptFakeTr);
 500       fOutputTrack->Add(hdistrFakeTr);
 501       fOutputTrack->Add(hd0f);
 502
 503     }
 504   }
 505
 506
 507   if(fOnOff[2] && fCuts->GetUseCentrality()){
 508
 509     //Centrality (Counters)
 510     fOutputCounters=new TList();
 511     fOutputCounters->SetOwner();
 512     fOutputCounters->SetName(GetOutputSlot(5)->GetContainer()->GetName());
 513
 514     AliCounterCollection *stdEstimator=new AliCounterCollection("stdEstimator");
 515     stdEstimator->AddRubric("run",500000);
 516     stdEstimator->AddRubric("centralityclass","-10_0/0_10/10_20/20_30/30_40/40_50/50_60/60_70/70_80/80_90/90_100/-990_-980");
 517     stdEstimator->Init();
 518     AliCounterCollection *secondEstimator=new AliCounterCollection("secondEstimator");
 519     secondEstimator->AddRubric("run",500000);
 520     secondEstimator->AddRubric("centralityclass","-10_0/0_10/10_20/20_30/30_40/40_50/50_60/60_70/70_80/80_90/90_100/-990_-980");
 521     secondEstimator->Init();
 522
 523     fOutputCounters->Add(stdEstimator);
 524     fOutputCounters->Add(secondEstimator);
 525
 526     //Centrality (Checks)
 527     fOutputCheckCentrality=new TList();
 528     fOutputCheckCentrality->SetOwner();
 529     fOutputCheckCentrality->SetName(GetOutputSlot(6)->GetContainer()->GetName());
 530
 531     TString hname="hNtrackletsIn";
 532     TH1F* hNtrackletsIn=new TH1F(hname.Data(),"Number of tracklets in Centrality range;ntracklets;Entries",5000,-0.5,4999.5);
 533
 534     hname="hMultIn";
 535     TH1F* hMultIn=new TH1F(hname.Data(),"Multiplicity;multiplicity in Centrality range;Entries",10000,-0.5,9999.5);
 536
 537     hname="hNtrackletsOut";
 538     TH1F* hNtrackletsOut=new TH1F(hname.Data(),"Number of tracklets out of Centrality range;ntracklets;Entries",5000,-0.5,4999.5);
 539
 540     hname="hMultOut";
 541     TH1F* hMultOut=new TH1F(hname.Data(),"Multiplicity out of Centrality range;multiplicity;Entries",10000,-0.5,9999.5);
 542
 543     hname="hMultvsPercentile";
 544     TH2F* hMultvsPercentile=new TH2F(hname.Data(),"Multiplicity vs Percentile;multiplicity;percentile",10000,-0.5,9999.5,240,-10.,110);
 545
 546     hname="hntrklvsPercentile";
 547     TH2F* hntrklvsPercentile=new TH2F(hname.Data(),"N tracklets vs Percentile;ntracklets;percentile",5000,-0.5,4999.5,240,-10.,110);
 548
 549     hname="hnTPCTracksvsPercentile";
 550     TH2F* hnTPCTracksvsPercentile=new TH2F(hname.Data(),"N TPC tracks vs Percentile;nTPCTracks;percentile",5000,-0.5,9999.5,240,-10.,110);
 551
 552     hname="hnTPCITSTracksvsPercentile";
 553     TH2F* hnTPCITSTracksvsPercentile=new TH2F(hname.Data(),"N TPC+ITS tracks vs Percentile;nTPCITSTracks;percentile",5000,-0.5,9999.5,240,-10.,110);
 554
 555     hname="hnTPCITS1SPDTracksvsPercentile";
 556     TH2F* hnTPCITS1SPDTracksvsPercentile=new TH2F(hname.Data(),"N TPC+ITS+1SPD tracks vs Percentile;nTPCITS1SPDTracks;percentile",5000,-0.5,9999.5,240,-10.,110);
 557
 558     hname="hV0MultiplicityPercentile";
 559     TH2F*hV0MultiplicityPercentile = new TH2F(hname.Data(),"V0 Multiplicity vs Percentile;V0 multiplicity;percentile",1000,-0.5,9999.5,120,-10.,110);
 560
 561     hname="hV0MultiplicityNtrackletsIn";
 562     TH2F*hV0MultiplicityNtrackletsIn = new TH2F(hname.Data(),"V0 Multiplicity vs Number of tracklets in the CC;V0 multiplicity;percentile",1000,-0.5,9999.5,5000,-0.5,4999.5);
 563
 564     hname="hStdPercentileSPDPercentile";
 565     TH2F* hStdPercentileSPDPercentile = new TH2F(hname.Data(),"Std estimator Percentile Vs SPD Percentile;Std estimator percentile;SPD percentile",120,-10.,110,120,-10.,110);
 566
 567     fOutputCheckCentrality->Add(hNtrackletsIn);
 568     fOutputCheckCentrality->Add(hNtrackletsOut);
 569     fOutputCheckCentrality->Add(hMultIn);
 570     fOutputCheckCentrality->Add(hMultOut);
 571     fOutputCheckCentrality->Add(hMultvsPercentile);
 572     fOutputCheckCentrality->Add(hntrklvsPercentile);
 573     fOutputCheckCentrality->Add(hnTPCTracksvsPercentile);
 574     fOutputCheckCentrality->Add(hnTPCITSTracksvsPercentile);
 575     fOutputCheckCentrality->Add(hnTPCITS1SPDTracksvsPercentile);
 576     fOutputCheckCentrality->Add(hV0MultiplicityPercentile);
 577     fOutputCheckCentrality->Add(hV0MultiplicityNtrackletsIn);
 578     fOutputCheckCentrality->Add(hStdPercentileSPDPercentile);
 579
 580     PostData(6,fOutputCheckCentrality);
 581
 582   } else{
 583     if(fOnOff[0]){
 584       TString hname="hNtracklets";
 585       TH1F* hNtracklets=new TH1F(hname.Data(),"Number of tracklets;ntracklets;Entries",5000,-0.5,4999.5);
 586
 587       hname="hMult";
 588       TH1F* hMult=new TH1F(hname.Data(),"Multiplicity;multiplicity;Entries",10000,-0.5,9999.5);
 589       fOutputTrack->Add(hNtracklets);
 590       fOutputTrack->Add(hMult);
 591     }
 592   }
 593
 594   //event selection (z vertex for the moment)
 595   if(fOnOff[3]){
 596     fOutputEvSelection=new TList();
 597     fOutputEvSelection->SetOwner();
 598     fOutputEvSelection->SetName(GetOutputSlot(7)->GetContainer()->GetName());
 599     AliCounterCollection *evselection=new AliCounterCollection("evselection");
 600     evselection->AddRubric("run",500000);
 601     evselection->AddRubric("evnonsel","zvtx");
 602     evselection->Init();
 603
 604     TH1F* hzvtx=new TH1F("hzvtx", "Distribution of z_{VTX};z_{VTX} [cm];Entries",100,-20,20);
 605
 606
 607     TH2F* hTrigCent=new TH2F("hTrigCent","Centrality vs. Trigger types",12,-1.5,10.5,12,-10,110);
 608     hTrigCent->GetXaxis()->SetBinLabel(1,"All");
 609     hTrigCent->GetXaxis()->SetBinLabel(2,"kAny");
 610     hTrigCent->GetXaxis()->SetBinLabel(3,"kMB");
 611     hTrigCent->GetXaxis()->SetBinLabel(4,"kINT7");
 612     hTrigCent->GetXaxis()->SetBinLabel(5,"kCINT5");
 613     hTrigCent->GetXaxis()->SetBinLabel(6,"kCent");
 614     hTrigCent->GetXaxis()->SetBinLabel(7,"kSemiCent");
 615     hTrigCent->GetXaxis()->SetBinLabel(8,"kEMC1+7");
 616     hTrigCent->GetXaxis()->SetBinLabel(9,"kEMCJET+GAMMA");
 617     hTrigCent->GetXaxis()->SetBinLabel(10,"Muons");
 618     hTrigCent->GetXaxis()->SetBinLabel(11,"PHOS");
 619     hTrigCent->GetXaxis()->SetBinLabel(12,"Others");
 620
 621     TH2F* hTrigMul=new TH2F("hTrigMul","Multiplicity vs. Trigger types",12,-1.5,10.5,100,0.,10000.);
 622     hTrigMul->GetXaxis()->SetBinLabel(1,"All");
 623     hTrigMul->GetXaxis()->SetBinLabel(2,"kAny");
 624     hTrigMul->GetXaxis()->SetBinLabel(3,"kMB");
 625     hTrigMul->GetXaxis()->SetBinLabel(4,"kINT7");
 626     hTrigMul->GetXaxis()->SetBinLabel(5,"kCINT5");
 627     hTrigMul->GetXaxis()->SetBinLabel(6,"kCent");
 628     hTrigMul->GetXaxis()->SetBinLabel(7,"kSemiCent");
 629     hTrigMul->GetXaxis()->SetBinLabel(8,"kEMC1+7");
 630     hTrigMul->GetXaxis()->SetBinLabel(9,"kEMCJET+GAMMA");
 631     hTrigMul->GetXaxis()->SetBinLabel(10,"Muons");
 632     hTrigMul->GetXaxis()->SetBinLabel(11,"PHOS");
 633     hTrigMul->GetXaxis()->SetBinLabel(12,"Others");
 634
 635     TH2F* hTrigCentSel=new TH2F("hTrigCentSel","Trigger types",12,-1.5,10.5,12,-10,110);
 636     hTrigCentSel->GetXaxis()->SetBinLabel(1,"All");
 637     hTrigCentSel->GetXaxis()->SetBinLabel(2,"kAny");
 638     hTrigCentSel->GetXaxis()->SetBinLabel(3,"kMB");
 639     hTrigCentSel->GetXaxis()->SetBinLabel(4,"kINT7");
 640     hTrigCentSel->GetXaxis()->SetBinLabel(5,"kCINT5");
 641     hTrigCentSel->GetXaxis()->SetBinLabel(6,"kCent");
 642     hTrigCentSel->GetXaxis()->SetBinLabel(7,"kSemiCent");
 643     hTrigCentSel->GetXaxis()->SetBinLabel(8,"kEMC1+7");
 644     hTrigCentSel->GetXaxis()->SetBinLabel(9,"kEMCJET+GAMMA");
 645     hTrigCentSel->GetXaxis()->SetBinLabel(10,"Muons");
 646     hTrigCentSel->GetXaxis()->SetBinLabel(11,"PHOS");
 647     hTrigCentSel->GetXaxis()->SetBinLabel(12,"Others");
 648
 649     AliCounterCollection *trigCounter=new AliCounterCollection("trigCounter");
 650     trigCounter->AddRubric("run",500000);
 651     trigCounter->AddRubric("triggerType","All/Any/MB/Cent/SemiCent/EMCAL/MUON/NoPhysSelMUON/NoPhysSelEvNot7/NoPhysSelCMUP1/NoPhysSelMB/NoPhysSelCent/NoPhysSelSemiCent");
 652     trigCounter->Init();
 653
 654     fOutputEvSelection->Add(evselection);
 655     fOutputEvSelection->Add(hzvtx);
 656     fOutputEvSelection->Add(hTrigCent);
 657     fOutputEvSelection->Add(hTrigMul);
 658     fOutputEvSelection->Add(hTrigCentSel);
 659     fOutputEvSelection->Add(trigCounter);
 660   }
 661   if(fOnOff[4]){ // FLOW OBSERVABLES
 662     fOutputFlowObs=new TList();
 663     fOutputFlowObs->SetOwner();
 664     fOutputFlowObs->SetName(GetOutputSlot(8)->GetContainer()->GetName());
 665
 666     fFlowEvent = new AliFlowEvent(3000);
 667     fRFPcuts = new AliFlowTrackCuts("rfpCuts");
 668
 669     TH2F *hFEvents = new TH2F("hFlowEvents","FlowEvent Selection",7,0,7,7,-10,60);
 670     hFEvents->GetXaxis()->SetBinLabel(1,"REACHED");
 671     hFEvents->GetXaxis()->SetBinLabel(2,"TRIGGERED");
 672     hFEvents->GetXaxis()->SetBinLabel(3,"kMB");
 673     hFEvents->GetXaxis()->SetBinLabel(4,"kCent");
 674     hFEvents->GetXaxis()->SetBinLabel(5,"kSemiC");
 675     hFEvents->GetXaxis()->SetBinLabel(6,"Triggered + vtx cut");
 676     hFEvents->GetXaxis()->SetBinLabel(7,"UnexpectedBehaviour");
 677     fOutputFlowObs->Add(hFEvents);
 678
 679     TProfile2D *hQ[3];
 680     TH2F *hAngleQ[3];
 681     TH3F *hPhiEta[3];
 682     TString ref[3] = {"FB1","FB128","VZE"};
 683     Int_t etabin[3] = {40,40,20};
 684     Int_t etamax[3] = { 1, 1, 5};
 685     for(Int_t i=0; i<3; ++i) {
 686       hQ[i]= new TProfile2D( Form("h%s_Q",ref[i].Data()),
 687                              Form("Q_{2} components for %s",ref[i].Data()),
 688                              4,0,4,12,0,60,"s");
 689       hQ[i]->GetXaxis()->SetBinLabel(1,"Qx^{-}");
 690       hQ[i]->GetXaxis()->SetBinLabel(2,"Qy^{-}");
 691       hQ[i]->GetXaxis()->SetBinLabel(3,"Qx^{+}");
 692       hQ[i]->GetXaxis()->SetBinLabel(4,"Qy^{+}");
 693       hQ[i]->GetYaxis()->SetTitle("Centrality");
 694       fOutputFlowObs->Add(hQ[i]);
 695
 696       hAngleQ[i] = new TH2F( Form("h%s_AngleQ",ref[i].Data()),
 697                              Form("#Psi_{2} for %s",ref[i].Data()),
 698                              72,0,TMath::Pi(),12,0,60);
 699       hAngleQ[i]->GetXaxis()->SetTitle( Form("#Psi_{2}^{%s}",ref[i].Data()) );
 700       hAngleQ[i]->GetYaxis()->SetTitle("Centrality");
 701       fOutputFlowObs->Add(hAngleQ[i]);
 702
 703       hPhiEta[i] = new TH3F( Form("h%s_PhiEta",ref[i].Data()),
 704                              Form("Eta vs Phi for %s",ref[i].Data()),
 705                              144,0,TMath::TwoPi(),etabin[i],-1.0*etamax[i],+1.0*etamax[i],12,0,60);
 706       hPhiEta[i]->GetXaxis()->SetTitle("Phi");
 707       hPhiEta[i]->GetYaxis()->SetTitle("Eta");
 708       hPhiEta[i]->GetZaxis()->SetTitle("Centrality");
 709       fOutputFlowObs->Add(hPhiEta[i]);
 710
 711     }
 712     TH3F *hTPCVZE_AngleQ = new TH3F("hTPCVZE_AngleQ","#Psi_{2}^{VZERO} vs #Psi_{2}^{TPC}",   72,0,TMath::Pi(),72,0,TMath::Pi(),12,0,60);
 713     hTPCVZE_AngleQ->GetXaxis()->SetTitle("#Psi_{2}^{TPC}");
 714     hTPCVZE_AngleQ->GetYaxis()->SetTitle("#Psi_{2}^{VZE}");
 715     hTPCVZE_AngleQ->GetZaxis()->SetTitle("Centrality");
 716     fOutputFlowObs->Add(hTPCVZE_AngleQ);
 717
 718     TH2F *hCentVsMultRPS = new TH2F("hCentVsMultRPS", " Centrality Vs. Multiplicity RPs",5000, 0, 5000.,12,0,60 );
 719     hCentVsMultRPS->GetXaxis()->SetTitle("Multiplicity RPs");
 720     hCentVsMultRPS->GetYaxis()->SetTitle("Centrality");
 721     fOutputFlowObs->Add(hCentVsMultRPS);
 722   }
 723 //  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 724 //  AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 725 //  AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 726 //  fCuts->GetPidHF()->SetPidResponse(pidResp);
 727   // Post the data
 728   PostData(1,fNEntries);
 729
 730   if(fOnOff[1]) PostData(2,fOutputPID);
 731   if(fOnOff[0]) PostData(3,fOutputTrack);
 732   PostData(4,fCuts);
 733   if(fOnOff[2]) PostData(5,fOutputCounters);
 734   if(fOnOff[3]) PostData(7,fOutputEvSelection);
 735   if(fOnOff[4]) PostData(8,fOutputFlowObs);
 736
 737   if(!fOnOff[0] && !fOnOff[1] && !fOnOff[2]) AliError("Nothing will be filled!");
 738 }
 739
 740 //___________________________________________________________________________
 741 void AliAnalysisTaskSEHFQA::UserExec(Option_t */*option*/)
 742 {
 743   // Execute analysis for current event
 744
 745   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 746   if(fDebug>2) printf("Analysing decay %d\n",fDecayChannel);
 747   // Post the data already here
 748   PostData(1,fNEntries);
 749   if(fOnOff[1]) PostData(2,fOutputPID);
 750   if(fOnOff[0]) PostData(3,fOutputTrack);
 751   PostData(4,fCuts);
 752   if(fOnOff[2]) {
 753     PostData(5,fOutputCounters);
 754     if(fCuts->GetUseCentrality()) PostData(6,fOutputCheckCentrality);
 755   }
 756
 757   TClonesArray *arrayProng =0;
 758   Int_t pdg=0;
 759   Int_t *pdgdaughters=0x0;
 760
 761   if(!aod && AODEvent() && IsStandardAOD()) {
 762     // In case there is an AOD handler writing a standard AOD, use the AOD
 763     // event in memory rather than the input (ESD) event.
 764     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 765     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 766     // have to taken from the AOD event hold by the AliAODExtension
 767     AliAODHandler* aodHandler = (AliAODHandler*)
 768       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 769     if(aodHandler->GetExtensions()) {
 770
 771       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 772       AliAODEvent *aodFromExt = ext->GetAOD();
 773
 774
 775
 776       switch(fDecayChannel){
 777       case 0:
 778         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 779         pdg=411;
 780         if(fReadMC){
 781           pdgdaughters =new Int_t[3];
 782           pdgdaughters[0]=211;//pi
 783           pdgdaughters[1]=321;//K
 784           pdgdaughters[2]=211;//pi
 785         }
 786         break;
 787       case 1:
 788         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 789         pdg=421;
 790         if(fReadMC){
 791           pdgdaughters =new Int_t[2];
 792           pdgdaughters[0]=211;//pi
 793           pdgdaughters[1]=321;//K
 794         }
 795         break;
 796       case 2:
 797         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 798         pdg=413;
 799         if(fReadMC){
 800           pdgdaughters =new Int_t[3];
 801           pdgdaughters[1]=211;//pi
 802           pdgdaughters[0]=321;//K
 803           pdgdaughters[2]=211;//pi (soft?)
 804         }
 805         break;
 806       case 3:
 807         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 808         pdg=431;
 809         if(fReadMC){
 810           pdgdaughters =new Int_t[3];
 811           pdgdaughters[0]=321;//K
 812           pdgdaughters[1]=321;//K
 813           pdgdaughters[2]=211;//pi
 814         }
 815         break;
 816       case 4:
 817         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 818         pdg=421;
 819         if(fReadMC){
 820           pdgdaughters =new Int_t[4];
 821           pdgdaughters[0]=321;
 822           pdgdaughters[1]=211;
 823           pdgdaughters[2]=211;
 824           pdgdaughters[3]=211;
 825         }
 826         break;
 827       case 5:
 828         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 829         pdg=4122;
 830         if(fReadMC){
 831           pdgdaughters =new Int_t[3];
 832           pdgdaughters[0]=2212;//p
 833           pdgdaughters[1]=321;//K
 834           pdgdaughters[2]=211;//pi
 835         }
 836         break;
 837       case kLambdactoV0:
 838         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadeHF");
 839         pdg=4122;
 840         if(fReadMC){
 841           pdgdaughters =new Int_t[3];
 842           pdgdaughters[0]=2212;//p
 843           pdgdaughters[1]=211;//pi
 844           pdgdaughters[2]=211;//pi
 845         }
 846         break;
 847       }
 848     }
 849   } else if(aod) {
 850     switch(fDecayChannel){
 851     case 0:
 852       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 853       pdg=411;
 854       if(fReadMC){
 855         pdgdaughters =new Int_t[3];
 856         pdgdaughters[0]=211;//pi
 857         pdgdaughters[1]=321;//K
 858         pdgdaughters[2]=211;//pi
 859       }
 860       break;
 861     case 1:
 862       arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 863       pdg=421;
 864       if(fReadMC){
 865         pdgdaughters =new Int_t[2];
 866         pdgdaughters[0]=211;//pi
 867         pdgdaughters[1]=321;//K
 868       }
 869       break;
 870     case 2:
 871       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
 872       pdg=413;
 873       if(fReadMC){
 874         pdgdaughters =new Int_t[3];
 875         pdgdaughters[1]=211;//pi
 876         pdgdaughters[0]=321;//K
 877         pdgdaughters[2]=211;//pi (soft?)
 878       }
 879       break;
 880     case 3:
 881       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 882       pdg=431;
 883       if(fReadMC){
 884         pdgdaughters =new Int_t[3];
 885         pdgdaughters[0]=321;//K
 886         pdgdaughters[1]=321;//K
 887         pdgdaughters[2]=211;//pi
 888       }
 889       break;
 890     case 4:
 891       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm4Prong");
 892       pdg=421;
 893       if(fReadMC){
 894         pdgdaughters =new Int_t[4];
 895         pdgdaughters[0]=321;
 896         pdgdaughters[1]=211;
 897         pdgdaughters[2]=211;
 898         pdgdaughters[3]=211;
 899       }
 900       break;
 901     case 5:
 902       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 903       pdg=4122;
 904       if(fReadMC){
 905         pdgdaughters =new Int_t[3];
 906         pdgdaughters[0]=2212;//p
 907         pdgdaughters[1]=321;//K
 908         pdgdaughters[2]=211;//pi
 909       }
 910       break;
 911       case kLambdactoV0:
 912         arrayProng=(TClonesArray*)aod->GetList()->FindObject("CascadeHF");
 913         pdg=4122;
 914         if(fReadMC){
 915           pdgdaughters =new Int_t[3];
 916           pdgdaughters[0]=2212;//p
 917           pdgdaughters[1]=211;//pi
 918           pdgdaughters[2]=211;//pi
 919         }
 920         break;
 921     }
 922   }
 923   Bool_t isSimpleMode=fSimpleMode;
 924   if(!arrayProng) {
 925     AliInfo("Branch not found! The output will contain only trak related histograms\n");
 926     isSimpleMode=kTRUE;
 927     fNEntries->Fill(2);
 928   }
 929
 930   TClonesArray *mcArray = 0;
 931   AliAODMCHeader *mcHeader = 0;
 932
 933   if(!aod) {
 934     delete [] pdgdaughters;
 935     return;
 936   }
 937
 938   //check if MC
 939   if(fReadMC) {
 940     // load MC particles
 941     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 942     if(!mcArray) {
 943       printf("AliAnalysisTaskSEHFQA::UserExec: MC particles branch not found!\n");
 944       delete [] pdgdaughters;
 945       return;
 946     }
 947
 948     // load MC header
 949     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 950     if(!mcHeader) {
 951       printf("AliAnalysisTaskSEHFQA::UserExec: MC header branch not found!\n");
 952       delete [] pdgdaughters;
 953       return;
 954     }
 955   }
 956
 957
 958   UInt_t evSelMask=((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 959   Double_t centrality=fCuts->GetCentrality(aod);
 960   Double_t multiplicity=aod->GetHeader()->GetRefMultiplicity();
 961   Int_t runNumber = aod->GetRunNumber();
 962   TString trigClass=aod->GetFiredTriggerClasses();
 963   Int_t nAODtracks=aod->GetNTracks();
 964   Int_t nSelTracksTPCOnly=0;
 965   Int_t nSelTracksTPCITS=0;
 966   Int_t nSelTracksTPCITS1SPD=0;
 967
 968   for (Int_t k=0;k<nAODtracks;k++){
 969     AliAODTrack* track=aod->GetTrack(k);
 970     if(track->GetID()<0) continue;
 971     Int_t nclsTot=0,nclsSPD=0;
 972     for(Int_t l=0;l<6;l++) {
 973       if(TESTBIT(track->GetITSClusterMap(),l)) {
 974         nclsTot++; if(l<2) nclsSPD++;
 975       }
 976     }
 977     UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
 978     if(TMath::Abs(track->Eta())<0.8 && nTPCClus>=70 && track->GetStatus()&AliESDtrack::kTPCrefit){
 979       if(track->TestFilterBit(1))  nSelTracksTPCOnly++;
 980       if(track->GetStatus()&AliESDtrack::kITSrefit){
 981         nSelTracksTPCITS++;
 982         if(nclsSPD>0) nSelTracksTPCITS1SPD++;
 983       }
 984     }
 985   }
 986
 987   if(fOnOff[4]) {
 988     FillFlowObs(aod);
 989     PostData(8,fOutputFlowObs);
 990   }
 991   if(fOnOff[3]){
 992     TH2F* hTrigC=(TH2F*)fOutputEvSelection->FindObject("hTrigCent");
 993     TH2F* hTrigM=(TH2F*)fOutputEvSelection->FindObject("hTrigMul");
 994     AliCounterCollection* trigCount=(AliCounterCollection*)fOutputEvSelection->FindObject("trigCounter");
 995
 996     hTrigC->Fill(-1.,centrality);
 997     hTrigM->Fill(-1.,multiplicity);
 998     trigCount->Count(Form("triggerType:All/Run:%d",runNumber));
 999     if(evSelMask==0){
1000       if(aod->GetEventType()!=7){
1001         trigCount->Count(Form("triggerType:NoPhysSelEvNot7/Run:%d",runNumber));
1002       }else if(trigClass.Contains("CMUP1")){
1003         trigCount->Count(Form("triggerType:NoPhysSelCMUP1/Run:%d",runNumber));
1004       }else if(trigClass.Contains("MUON")){
1005         trigCount->Count(Form("triggerType:NoPhysSelMUON/Run:%d",runNumber));
1006       }else if(trigClass.Contains("CPBI2_B1-B") || trigClass.Contains(" CPBI2WU_B1-B")){
1007         trigCount->Count(Form("triggerType:NoPhysSelMB/Run:%d",runNumber));
1008       }else if(trigClass.Contains("CCENT") || trigClass.Contains("CVHN")){
1009         trigCount->Count(Form("triggerType:NoPhysSelCent/Run:%d",runNumber));
1010       }else if(trigClass.Contains("CSEMI") || trigClass.Contains("CVLN")){
1011         trigCount->Count(Form("triggerType:NoPhysSelSemiCent/Run:%d",runNumber));
1012       }
1013     }
1014     if(evSelMask & AliVEvent::kAny){
1015       hTrigC->Fill(0.,centrality);
1016       hTrigM->Fill(0.,multiplicity);
1017       trigCount->Count(Form("triggerType:Any/Run:%d",runNumber));
1018     }
1019     if(evSelMask & AliVEvent::kMB){
1020       hTrigC->Fill(1.,centrality);
1021       hTrigM->Fill(1.,multiplicity);
1022       trigCount->Count(Form("triggerType:MB/Run:%d",runNumber));
1023     }
1024     if(evSelMask & AliVEvent::kINT7){
1025       hTrigC->Fill(2.,centrality);
1026       hTrigM->Fill(2.,multiplicity);
1027     }
1028     if(evSelMask & AliVEvent::kCINT5){
1029       hTrigC->Fill(3.,centrality);
1030       hTrigM->Fill(3.,multiplicity);
1031     }
1032     if(evSelMask & AliVEvent::kCentral){
1033       hTrigC->Fill(4.,centrality);
1034       hTrigM->Fill(4.,multiplicity);
1035       trigCount->Count(Form("triggerType:Cent/Run:%d",runNumber));
1036     }
1037     if(evSelMask & AliVEvent::kSemiCentral){
1038       hTrigC->Fill(5.,centrality);
1039       hTrigM->Fill(5.,multiplicity);
1040       trigCount->Count(Form("triggerType:SemiCent/Run:%d",runNumber));
1041     }
1042     if(evSelMask & (AliVEvent::kEMC1 | AliVEvent::kEMC7)){
1043       hTrigC->Fill(6.,centrality);
1044       hTrigM->Fill(6.,multiplicity);
1045     }
1046     if(evSelMask & (AliVEvent::kEMCEJE | AliVEvent::kEMCEGA)){
1047       hTrigC->Fill(7.,centrality);
1048       hTrigM->Fill(7.,multiplicity);
1049       trigCount->Count(Form("triggerType:EMCAL/Run:%d",runNumber));
1050     }
1051     if(evSelMask & (((AliVEvent::kCMUS5 | AliVEvent::kMUSH7) | (AliVEvent::kMUL7 | AliVEvent::kMUU7)) |  (AliVEvent::kMUS7 | AliVEvent::kMUON))){
1052       hTrigC->Fill(8.,centrality);
1053       hTrigM->Fill(8.,multiplicity);
1054       trigCount->Count(Form("triggerType:MUON/Run:%d",runNumber));
1055     }
1056     if(evSelMask & (AliVEvent::kPHI1 | AliVEvent::kPHI7)){
1057       hTrigC->Fill(9.,centrality);
1058       hTrigM->Fill(9.,multiplicity);
1059     }
1060     if(evSelMask & (AliVEvent::kDG5 | AliVEvent::kZED)){
1061       hTrigC->Fill(10.,centrality);
1062       hTrigM->Fill(10.,multiplicity);
1063     }
1064   }
1065
1066
1067   // fix for temporary bug in ESDfilter
1068   // the AODs with null vertex pointer didn't pass the PhysSel
1069   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) {
1070     delete [] pdgdaughters;
1071     return;
1072   }
1073
1074   // count event
1075   fNEntries->Fill(0);
1076
1077   //count events with good vertex
1078   // AOD primary vertex
1079   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
1080
1081   Double_t pos[3],cov[6];
1082   vtx1->GetXYZ(pos);
1083   vtx1->GetCovarianceMatrix(cov);
1084   const AliESDVertex vESD(pos,cov,100.,100);
1085
1086   TString primTitle = vtx1->GetTitle();
1087   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) fNEntries->Fill(4);
1088
1089   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
1090   //TString trigclass=aod->GetFiredTriggerClasses();
1091   //if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNEntries->Fill(5); //tmp
1092
1093
1094
1095
1096   Bool_t evSelbyCentrality=kTRUE,evSelected=kTRUE,evSelByVertex=kTRUE,evselByPileup=kTRUE,evSelByPS=kTRUE;
1097   //select event
1098   if(!fCuts->IsEventSelected(aod)) {
1099     evSelected=kFALSE;
1100     if(fCuts->GetWhyRejection()==1) {fNEntries->Fill(1); evselByPileup=kFALSE;}// rejected for pileup
1101     if(fCuts->GetWhyRejection()==2 || fCuts->GetWhyRejection()==3) evSelbyCentrality=kFALSE; //rejected by centrality
1102     if(fCuts->GetWhyRejection()==4) evSelByVertex=kFALSE; //rejected by vertex
1103     if(fCuts->GetWhyRejection()==5) fNEntries->Fill(5);//tmp
1104     if(fCuts->GetWhyRejection()==6 && fOnOff[3]) ((AliCounterCollection*)fOutputEvSelection->FindObject("evselection"))->Count(Form("evnonsel:zvtx/Run:%d",runNumber));
1105     if(fCuts->GetWhyRejection()==7) { evSelByPS=kFALSE; }
1106   }
1107   if(evSelected){
1108     TH2F* hTrigS=(TH2F*)fOutputEvSelection->FindObject("hTrigCentSel");
1109     hTrigS->Fill(-1.,centrality);
1110
1111     if(evSelMask & AliVEvent::kAny) hTrigS->Fill(0.,centrality);
1112     if(evSelMask & AliVEvent::kMB) hTrigS->Fill(1.,centrality);
1113     if(evSelMask & AliVEvent::kINT7) hTrigS->Fill(2.,centrality);
1114     if(evSelMask & AliVEvent::kCINT5) hTrigS->Fill(3.,centrality);
1115     if(evSelMask & AliVEvent::kCentral) hTrigS->Fill(4.,centrality);
1116     if(evSelMask & AliVEvent::kSemiCentral) hTrigS->Fill(5.,centrality);
1117     if(evSelMask & (AliVEvent::kEMC1 | AliVEvent::kEMC7)) hTrigS->Fill(6.,centrality);
1118     if(evSelMask & (AliVEvent::kEMCEJE | AliVEvent::kEMCEGA)) hTrigS->Fill(7.,centrality);
1119     if(evSelMask & (((AliVEvent::kCMUS5 | AliVEvent::kMUSH7) | (AliVEvent::kMUL7 | AliVEvent::kMUU7)) |  (AliVEvent::kMUS7 | AliVEvent::kMUON))) hTrigS->Fill(8.,centrality);
1120     if(evSelMask & (AliVEvent::kPHI1 | AliVEvent::kPHI7)) hTrigS->Fill(9.,centrality);
1121     if(evSelMask & (AliVEvent::kDG5 | AliVEvent::kZED)) hTrigS->Fill(10.,centrality);
1122   }
1123
1124   if(evSelected || (!evSelbyCentrality && evSelByVertex && evselByPileup && evSelByPS)){ //events selected or not selected because of centrality
1125     if(fOnOff[2] && fCuts->GetUseCentrality()){
1126
1127       Float_t stdCentf=fCuts->GetCentrality(aod);
1128       Int_t stdCent = (Int_t)(stdCentf+0.5);
1129       Float_t secondCentf =fCuts->GetCentrality(aod,fEstimator);
1130       Int_t secondCent = (Int_t)(secondCentf+0.5);
1131       Int_t mincent=stdCent-stdCent%10;
1132       AliAODVZERO *vzeroAOD = (AliAODVZERO*)aod->GetVZEROData();
1133       Float_t vzeroMult = vzeroAOD->GetMTotV0A() +  vzeroAOD->GetMTotV0C();
1134       AliCentrality *aodcent = aod->GetCentrality();
1135       Float_t spdCentf = aodcent->GetCentralityPercentile("CL1");
1136       if(stdCentf==-1) {
1137         mincent=-10;
1138         stdCent=-1;
1139       }
1140       if(mincent==100)mincent--;
1141       ((AliCounterCollection*)fOutputCounters->FindObject("stdEstimator"))->Count(Form("centralityclass:%d_%d/Run:%d",mincent,mincent+10,runNumber));
1142
1143       mincent=secondCent-secondCent%10;
1144       if(secondCentf==-1) {
1145         mincent=-10;
1146         secondCent=-1;
1147       }
1148       if(mincent==100)mincent--;
1149       ((AliCounterCollection*)fOutputCounters->FindObject("secondEstimator"))->Count(Form("centralityclass:%d_%d/Run:%d",mincent,mincent+10,runNumber));
1150
1151       if(stdCent<fCuts->GetMinCentrality() || stdCent>fCuts->GetMaxCentrality()){
1152         ((TH1F*)fOutputCheckCentrality->FindObject("hNtrackletsOut"))->Fill(aod->GetTracklets()->GetNumberOfTracklets());
1153         ((TH1F*)fOutputCheckCentrality->FindObject("hMultOut"))->Fill(aod->GetHeader()->GetRefMultiplicity());
1154       }else{
1155         ((TH1F*)fOutputCheckCentrality->FindObject("hNtrackletsIn"))->Fill(aod->GetTracklets()->GetNumberOfTracklets());
1156         ((TH1F*)fOutputCheckCentrality->FindObject("hMultIn"))->Fill(aod->GetHeader()->GetRefMultiplicity());
1157       }
1158       ((TH2F*)fOutputCheckCentrality->FindObject("hMultvsPercentile"))->Fill(aod->GetHeader()->GetRefMultiplicity(),stdCentf);
1159       ((TH2F*)fOutputCheckCentrality->FindObject("hntrklvsPercentile"))->Fill(aod->GetTracklets()->GetNumberOfTracklets(),stdCentf);
1160       ((TH2F*)fOutputCheckCentrality->FindObject("hnTPCTracksvsPercentile"))->Fill(nSelTracksTPCOnly,stdCentf);
1161       ((TH2F*)fOutputCheckCentrality->FindObject("hnTPCITSTracksvsPercentile"))->Fill(nSelTracksTPCITS,stdCentf);
1162       ((TH2F*)fOutputCheckCentrality->FindObject("hnTPCITS1SPDTracksvsPercentile"))->Fill(nSelTracksTPCITS1SPD,stdCentf);
1163       ((TH2F*)fOutputCheckCentrality->FindObject("hV0MultiplicityPercentile"))->Fill(vzeroMult,stdCentf);
1164       ((TH2F*)fOutputCheckCentrality->FindObject("hV0MultiplicityNtrackletsIn"))->Fill(vzeroMult,aod->GetTracklets()->GetNumberOfTracklets());
1165       ((TH2F*)fOutputCheckCentrality->FindObject("hStdPercentileSPDPercentile"))->Fill(stdCentf,spdCentf);
1166
1167       PostData(6,fOutputCheckCentrality);
1168
1169     } else{
1170       if(fOnOff[0]){
1171         ((TH1F*)fOutputTrack->FindObject("hNtracklets"))->Fill(aod->GetTracklets()->GetNumberOfTracklets());
1172         ((TH1F*)fOutputTrack->FindObject("hMult"))->Fill(aod->GetHeader()->GetRefMultiplicity());
1173       }
1174     }
1175   }
1176
1177   if(fOnOff[3]){
1178     const AliVVertex *vertex = aod->GetPrimaryVertex();
1179     Double_t zvtx=vertex->GetZ();
1180     if(evSelected || (!evSelected && TMath::Abs(zvtx) > 10.))
1181     ((TH1F*)fOutputEvSelection->FindObject("hzvtx"))->Fill(zvtx);
1182   }
1183
1184   if(!evSelected) {
1185     delete [] pdgdaughters;
1186     return; //discard all events not selected (vtx and/or centrality)
1187   }
1188
1189
1190   AliAODPidHF* pidHF=fCuts->GetPidHF();
1191   if(!pidHF) {
1192     delete [] pdgdaughters;
1193     return;
1194   }
1195   //AliPIDResponse* respF=pidHF->GetPidResponse();
1196   AliTPCPIDResponse* tpcres=new AliTPCPIDResponse();
1197   Bool_t oldPID=pidHF->GetOldPid();
1198   if(oldPID){
1199     Double_t alephParameters[5];
1200     pidHF->GetTPCBetheBlochParams(alephParameters);
1201     tpcres->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
1202   }
1203
1204
1205   Int_t ntracks=0;
1206   Int_t isGoodTrack=0, isFakeTrack=0, isSelTrack=0;
1207
1208   if(aod) ntracks=aod->GetNTracks();
1209
1210   if(fOnOff[0] || fOnOff[1]){
1211     //loop on tracks in the event
1212     for (Int_t k=0;k<ntracks;k++){
1213       AliAODTrack* track=aod->GetTrack(k);
1214       if(track->GetID()<0) continue;
1215       AliAODPid *pid = track->GetDetPid();
1216       if(!pid && fDebug>1) cout<<"No AliAODPid found"<<endl;
1217
1218       if(pid && fOnOff[1]){
1219         Double_t times[AliPID::kSPECIES];
1220         pid->GetIntegratedTimes(times);
1221
1222         Double_t tofRes[AliPID::kSPECIES];
1223         pid->GetTOFpidResolution(tofRes);
1224
1225         //check TOF
1226         TH1F* htmpfl=((TH1F*)fOutputPID->FindObject("hTOFflags"));
1227         htmpfl->Fill(0.);
1228         if (track->GetStatus()&AliESDtrack::kTPCout) htmpfl->Fill(1.);
1229         if (track->GetStatus()&AliESDtrack::kTOFout) htmpfl->Fill(2.);
1230         if (track->GetStatus()&AliESDtrack::kTIME) htmpfl->Fill(3.);
1231         if (track->GetStatus()&AliESDtrack::kTOFpid) htmpfl->Fill(4.);
1232         if (track->GetStatus()&AliESDtrack::kTOFmismatch) htmpfl->Fill(5.);
1233
1234         if(pidHF && pidHF->CheckStatus(track,"TOF")){
1235           ((TH1F*)fOutputPID->FindObject("hTOFtime"))->Fill(times[AliPID::kProton]);
1236           ((TH2F*)fOutputPID->FindObject("hTOFtimeKaonHyptime"))->Fill(track->P(),pid->GetTOFsignal()-times[3]); //3 is kaon
1237           ((TH1F*)fOutputPID->FindObject("hTOFsig"))->Fill(pid->GetTOFsignal());
1238           if (pid->GetTOFsignal()< 0) ((TH1F*)fOutputPID->FindObject("hTOFsig"))->Fill(-1);
1239
1240           Double_t nsigma[3]={-10,-10,-10};
1241           pidHF->GetnSigmaTOF(track,(Int_t)AliPID::kPion,nsigma[0]);
1242           pidHF->GetnSigmaTOF(track,(Int_t)AliPID::kKaon,nsigma[1]);
1243           pidHF->GetnSigmaTOF(track,(Int_t)AliPID::kProton,nsigma[2]);
1244
1245           ((TH2F*)fOutputPID->FindObject("hTOFsigmaKSigPid"))->Fill(track->P(),nsigma[1]);
1246           ((TH2F*)fOutputPID->FindObject("hTOFsigmaPionSigPid"))->Fill(track->P(),nsigma[0]);
1247           ((TH2F*)fOutputPID->FindObject("hTOFsigmaProtonSigPid"))->Fill(track->P(),nsigma[2]);
1248           if(fReadMC){
1249             Int_t label=track->GetLabel();
1250             if(label<=0) continue;
1251             AliMCParticle* mcpart=(AliMCParticle*)mcArray->At(label);
1252             if(mcpart){
1253               Int_t abspdgcode=TMath::Abs(mcpart->PdgCode());
1254               if(abspdgcode==211) ((TH2F*)fOutputPID->FindObject("hTOFsigmaMCPionSigPid"))->Fill(track->P(),nsigma[0]);
1255               if(abspdgcode==321) ((TH2F*)fOutputPID->FindObject("hTOFsigmaMCKSigPid"))->Fill(track->P(),nsigma[1]);
1256               if(abspdgcode==2212) ((TH2F*)fOutputPID->FindObject("hTOFsigmaMCProtonSigPid"))->Fill(track->P(),nsigma[2]);
1257
1258             }
1259           }
1260
1261           for (Int_t iS=2; iS<5; iS++){ //we plot TOF Pid resolution for 3-sigma identified particles
1262             if ( TMath::Abs(nsigma[iS-2])<3.){
1263               switch (iS) {
1264               case AliPID::kPion:
1265                 ((TH1F*)fOutputPID->FindObject("hTOFsigPid3sigPion"))->Fill(tofRes[iS]);
1266                 break;
1267               case AliPID::kKaon:
1268                 ((TH1F*)fOutputPID->FindObject("hTOFsigPid3sigKaon"))->Fill(tofRes[iS]);
1269                 break;
1270               case AliPID::kProton:
1271                 ((TH1F*)fOutputPID->FindObject("hTOFsigPid3sigProton"))->Fill(tofRes[iS]);
1272                 break;
1273               default:
1274                 break;
1275               }
1276             }
1277           }
1278         }//if TOF status
1279         //}
1280
1281         if(pidHF && pidHF->CheckStatus(track,"TPC")){
1282
1283           Double_t TPCp=pid->GetTPCmomentum();
1284           Double_t TPCsignal=pid->GetTPCsignal();
1285           ((TH1F*)fOutputPID->FindObject("hTPCsig"))->Fill(TPCsignal);
1286           ((TH1F*)fOutputPID->FindObject("hTPCsigvsp"))->Fill(TPCp,TPCsignal);
1287           //if (pidHF->IsKaonRaw(track, "TOF"))
1288           Double_t nsigma[3]={-10,-10,-10};
1289           pidHF->GetnSigmaTPC(track,(Int_t)AliPID::kPion,nsigma[0]);
1290           pidHF->GetnSigmaTPC(track,(Int_t)AliPID::kKaon,nsigma[1]);
1291           pidHF->GetnSigmaTPC(track,(Int_t)AliPID::kProton,nsigma[2]);
1292
1293           ((TH2F*)fOutputPID->FindObject("hTPCsigmaK"))->Fill(TPCp,nsigma[1]);
1294
1295           ((TH2F*)fOutputPID->FindObject("hTPCsigmaPion"))->Fill(TPCp,nsigma[0]);
1296           ((TH2F*)fOutputPID->FindObject("hTPCsigmaProton"))->Fill(TPCp,nsigma[2]);
1297
1298           if(fReadMC){
1299             Int_t label=track->GetLabel();
1300             if(label<=0) continue;
1301             AliMCParticle* mcpart=(AliMCParticle*)mcArray->At(label);
1302             if(mcpart){
1303               Int_t abspdgcode=TMath::Abs(mcpart->PdgCode());
1304               if(abspdgcode==211) ((TH2F*)fOutputPID->FindObject("hTPCsigmaMCPion"))->Fill(track->P(),nsigma[0]);
1305               if(abspdgcode==321) ((TH2F*)fOutputPID->FindObject("hTPCsigmaMCK"))->Fill(track->P(),nsigma[1]);
1306               if(abspdgcode==2212) ((TH2F*)fOutputPID->FindObject("hTPCsigmaMCProton"))->Fill(track->P(),nsigma[2]);
1307
1308             }
1309
1310           }
1311
1312
1313         }//if TPC status
1314       } //end PID histograms
1315
1316       Int_t nclsTot=0,nclsSPD=0;
1317
1318       //check clusters of the tracks
1319       if(fOnOff[0]){
1320
1321         ((TH1F*)fOutputTrack->FindObject("hnLayerITS"))->Fill(-1);
1322         for(Int_t l=0;l<6;l++) {
1323           if(TESTBIT(track->GetITSClusterMap(),l)) {
1324             ((TH1F*)fOutputTrack->FindObject("hnLayerITS"))->Fill(l);
1325             nclsTot++; if(l<2) nclsSPD++;
1326           }
1327         }
1328         ((TH1F*)fOutputTrack->FindObject("hnClsITS"))->Fill(nclsTot);
1329         ((TH1F*)fOutputTrack->FindObject("hnClsSPD"))->Fill(nclsSPD);
1330         if(track->Pt()>0.3 &&
1331            TMath::Abs(track->Eta())<0.8 &&
1332            track->GetStatus()&AliESDtrack::kITSrefit &&
1333            track->GetStatus()&AliESDtrack::kTPCrefit &&
1334            nclsSPD>0){
1335           ((TH1F*)fOutputTrack->FindObject("hnClsITSselTr"))->Fill(nclsTot);
1336         }
1337         if(!(track->GetStatus()&AliESDtrack::kTPCin) && track->GetStatus()&AliESDtrack::kITSrefit && !(track->GetStatus()&AliESDtrack::kITSpureSA)){//tracks retrieved in the ITS and not reconstructed in the TPC
1338           ((TH1F*)fOutputTrack->FindObject("hnClsITS-SA"))->Fill(nclsTot);
1339           ((TH1F*)fOutputTrack->FindObject("hnLayerITS"))->Fill(-1);
1340           for(Int_t l=0;l<6;l++) {
1341             if(TESTBIT(track->GetITSClusterMap(),l)) {
1342               ((TH1F*)fOutputTrack->FindObject("hnLayerITSsa"))->Fill(l);
1343             }
1344           }
1345         }
1346         Int_t label=0;
1347         if(fReadMC){
1348           label=track->GetLabel();
1349           if (label<0)fNEntries->Fill(8);
1350           else fNEntries->Fill(9);
1351         }
1352
1353
1354         if (track->Pt()>0.3 &&
1355             track->GetStatus()&AliESDtrack::kTPCrefit &&
1356             track->GetStatus()&AliESDtrack::kITSrefit &&
1357             /*nclsTot>3 &&*/
1358             nclsSPD>0) {//count good tracks
1359
1360
1361           if(fReadMC && label<0) {
1362             ((TH1F*)fOutputTrack->FindObject("hptFakeTr"))->Fill(track->Pt());
1363             isFakeTrack++;
1364           } else {
1365             ((TH1F*)fOutputTrack->FindObject("hptGoodTr"))->Fill(track->Pt());
1366             isGoodTrack++;
1367           }
1368         }
1369         if(fCuts->IsDaughterSelected(track,&vESD,fCuts->GetTrackCuts())){
1370           isSelTrack++;
1371         }//select tracks for our analyses
1372       } //fill track histos
1373     } //end loop on tracks
1374
1375     //fill once per event
1376     if(fOnOff[0]){
1377       if (fReadMC) ((TH1F*)fOutputTrack->FindObject("hdistrFakeTr"))->Fill(isFakeTrack);
1378       ((TH1F*)fOutputTrack->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
1379       ((TH1F*)fOutputTrack->FindObject("hdistrSelTr"))->Fill(isSelTrack);
1380     }
1381
1382     if(!isSimpleMode){
1383       // loop over candidates
1384       Int_t nCand = arrayProng->GetEntriesFast();
1385       Int_t ndaugh=3;
1386       if(fDecayChannel==AliAnalysisTaskSEHFQA::kD0toKpi) ndaugh=2;
1387       if(fDecayChannel==AliAnalysisTaskSEHFQA::kD0toKpipipi) ndaugh=4;
1388
1389       for (Int_t iCand = 0; iCand < nCand; iCand++) {
1390         AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayProng->UncheckedAt(iCand);
1391         if(d->GetSelectionMap()) {
1392           if(fDecayChannel==AliAnalysisTaskSEHFQA::kD0toKpi && !d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)) continue; //skip the D0 from Dstar
1393           if(fDecayChannel==AliAnalysisTaskSEHFQA::kDplustoKpipi && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)) continue; //skip the 3 prong !D+
1394         }
1395
1396         if(fReadMC){
1397           Int_t labD = d->MatchToMC(pdg,mcArray,ndaugh,pdgdaughters);
1398           if(labD>=0){
1399             AliAODMCParticle *partD = (AliAODMCParticle*)mcArray->At(labD);
1400             Int_t label=partD->GetMother();
1401             AliAODMCParticle *mot = (AliAODMCParticle*)mcArray->At(label);
1402             while(label>=0){//get first mother
1403               mot = (AliAODMCParticle*)mcArray->At(label);
1404               label=mot->GetMother();
1405             }
1406             if(mot){
1407               Int_t pdgMotCode = mot->GetPdgCode();
1408
1409               if(TMath::Abs(pdgMotCode)==4) fNEntries->Fill(6); //from primary charm
1410               if(TMath::Abs(pdgMotCode)==5) fNEntries->Fill(7); //from beauty
1411             }
1412           }
1413         }//end MC
1414         else fNEntries->Fill(6); //count the candidates (data)
1415
1416         for(Int_t id=0;id<ndaugh;id++){
1417
1418           //other histograms to be filled when the cut object is given
1419           AliAODTrack* track=(AliAODTrack*)d->GetDaughter(id);
1420
1421           //track quality
1422
1423           if (fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(pdg)) && fCuts->IsSelected(d,AliRDHFCuts::kTracks,aod)) {
1424
1425             Int_t label=0;
1426             if(fReadMC)label=track->GetLabel();
1427             if(fOnOff[0]){
1428
1429               if(fReadMC && label<0) {
1430                 isFakeTrack++;
1431                 ((TH1F*)fOutputTrack->FindObject("hptFakeTrFromDaugh"))->Fill(track->Pt());
1432
1433                 ((TH1F*)fOutputTrack->FindObject("hd0f"))->Fill(d->Getd0Prong(id));
1434               } else {
1435                 ((TH1F*)fOutputTrack->FindObject("hptGoodTrFromDaugh"))->Fill(track->Pt());
1436                 ((TH1F*)fOutputTrack->FindObject("hd0"))->Fill(d->Getd0Prong(id));
1437                 Double_t d0rphiz[2],covd0[3];
1438                 Bool_t isDCA=track->PropagateToDCA(aod->GetPrimaryVertex(),aod->GetMagneticField(),9999.,d0rphiz,covd0);
1439                 if(isDCA) ((TH1F*)fOutputTrack->FindObject("hd0z"))->Fill(d0rphiz[1]);
1440               }
1441             }
1442             if (fCuts->IsSelected(d,AliRDHFCuts::kAll,aod) && fOnOff[1]){
1443
1444               AliAODPid *pid = track->GetDetPid();
1445               if(pid){
1446                 Double_t times[5];
1447                 pid->GetIntegratedTimes(times);
1448                 if(pidHF && pidHF->CheckStatus(track,"TOF")) ((TH2F*)fOutputPID->FindObject("hTOFtimeKaonHyptimeAC"))->Fill(track->P(),pid->GetTOFsignal()-times[AliPID::kKaon]);
1449                 if(pidHF && pidHF->CheckStatus(track,"TPC")) ((TH2F*)fOutputPID->FindObject("hTPCsigvspAC"))->Fill(pid->GetTPCmomentum(),pid->GetTPCsignal());
1450               }
1451               fNEntries->Fill(3);
1452             } //end analysis cuts
1453           } //end acceptance and track cuts
1454         } //end loop on tracks in the candidate
1455       } //end loop on candidates
1456
1457     }
1458   } //end if on pid or track histograms
1459
1460   delete tpcres;
1461   delete [] pdgdaughters;
1462   PostData(1,fNEntries);
1463   if(fOnOff[1]) PostData(2,fOutputPID);
1464   if(fOnOff[0]) PostData(3,fOutputTrack);
1465   PostData(4,fCuts);
1466   if(fOnOff[2]) PostData(5,fOutputCounters);
1467   //Post data 6 done in case of centrality on
1468
1469 }
1470
1471 //____________________________________________________________________________
1472 void AliAnalysisTaskSEHFQA::FillFlowObs(AliAODEvent *aod){
1473   //fills the flow observables
1474   Double_t cc;
1475   cc = fCuts->GetCentrality(aod);
1476   ((TH2F*) fOutputFlowObs->FindObject("hFlowEvents"))->Fill(0., cc);
1477
1478   UInt_t mask=((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1479   UInt_t trigger=AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral;
1480   if(mask & trigger) {
1481     ((TH2F*) fOutputFlowObs->FindObject("hFlowEvents"))->Fill(1.,cc); // fired
1482     if (mask & AliVEvent::kMB) ((TH2F*) fOutputFlowObs->FindObject("hFlowEvents"))->Fill(2.,cc);
1483     if (mask & AliVEvent::kCentral) ((TH2F*) fOutputFlowObs->FindObject("hFlowEvents"))->Fill(3.,cc);
1484     if (mask & AliVEvent::kSemiCentral) ((TH2F*) fOutputFlowObs->FindObject("hFlowEvents"))->Fill(4.,cc);
1485     Bool_t rejected=false;
1486     if(cc<0 || cc>60) rejected=true;
1487     const AliVVertex *vertex = aod->GetPrimaryVertex();
1488     Double_t zvtx=vertex->GetZ();
1489     if(TMath::Abs(zvtx)>fCuts->GetMaxVtxZ()) rejected=true;
1490     if(rejected) return; //not interesting for flow QA
1491   } else {
1492     return;
1493   }
1494
1495   // event accepted
1496   ((TH2F*) fOutputFlowObs->FindObject("hFlowEvents"))->Fill(5.,cc);
1497   fRFPcuts->SetParamType(AliFlowTrackCuts::kGlobal);
1498   fRFPcuts->SetPtRange(0.2,5.);
1499   fRFPcuts->SetEtaRange(-0.8,0.8);
1500   fRFPcuts->SetMinNClustersTPC(70);
1501   fRFPcuts->SetMinChi2PerClusterTPC(0.2);
1502   fRFPcuts->SetMaxChi2PerClusterTPC(4.0);
1503   fRFPcuts->SetAcceptKinkDaughters(kFALSE);
1504   fRFPcuts->SetEvent(aod);
1505
1506   TString ref[3] = {"FB1","FB128","VZE"};
1507   Double_t psi[3];
1508   for(Int_t i=0; i!=3; ++i) {
1509     if(i==0) { // switching to bit 1
1510       fRFPcuts->SetMinimalTPCdedx(10.);
1511       fRFPcuts->SetAODfilterBit(1);
1512     } else { // switching to bit 128
1513       fRFPcuts->SetMinimalTPCdedx(-1);
1514       fRFPcuts->SetAODfilterBit(128);
1515     }
1516     if(i>1) {
1517       fRFPcuts->SetParamType(AliFlowTrackCuts::kV0);
1518       fRFPcuts->SetEtaRange(-5,+5);
1519       fRFPcuts->SetPhiMin(0);
1520       fRFPcuts->SetPhiMax(TMath::TwoPi());
1521     }
1522     fFlowEvent->Fill(fRFPcuts,fRFPcuts);
1523     fFlowEvent->TagSubeventsInEta(-5,0,0,+5);
1524     // getting informationt
1525     AliFlowVector vQ, vQaQb[2];
1526     fFlowEvent->Get2Qsub(vQaQb,2);
1527     vQ = vQaQb[0]+vQaQb[1];
1528     Double_t dMa=vQaQb[0].GetMult();
1529     Double_t dMb=vQaQb[1].GetMult();
1530     if( dMa<2 || dMb<2 ) {
1531       ((TH2F*) fOutputFlowObs->FindObject("hFlowEvents"))->Fill(6.,cc); //???
1532       continue;
1533     }
1534     psi[i] = vQ.Phi()/2;
1535     // publishing
1536     ((TProfile2D*) fOutputFlowObs->FindObject( Form("h%s_Q",ref[i].Data())))->Fill(0,cc,vQaQb[0].X()/dMa,dMa); // Qx-
1537     ((TProfile2D*) fOutputFlowObs->FindObject( Form("h%s_Q",ref[i].Data())))->Fill(1,cc,vQaQb[0].Y()/dMa,dMa); // Qy-
1538     ((TProfile2D*) fOutputFlowObs->FindObject( Form("h%s_Q",ref[i].Data())))->Fill(2,cc,vQaQb[1].X()/dMb,dMb); // Qx+
1539     ((TProfile2D*) fOutputFlowObs->FindObject( Form("h%s_Q",ref[i].Data())))->Fill(3,cc,vQaQb[1].Y()/dMb,dMb); // Qy+
1540     ((TH2F*) fOutputFlowObs->FindObject( Form("h%s_AngleQ",ref[i].Data()) ))->Fill(psi[i],cc); // Psi
1541     AliFlowTrackSimple *track;
1542     for(Int_t t=0; t!=fFlowEvent->NumberOfTracks(); ++t) {
1543       track = (AliFlowTrackSimple*) fFlowEvent->GetTrack(t);
1544       if(!track) continue;
1545       if(!track->InRPSelection()) continue;
1546       ((TH3F*) fOutputFlowObs->FindObject( Form("h%s_PhiEta",ref[i].Data()) ))->Fill(track->Phi(),track->Eta(),cc,track->Weight()); //PhiEta
1547     }
1548
1549   //histo filled only for TPCFB1
1550   if (i==0) {
1551     ((TH2F*) fOutputFlowObs->FindObject("hCentVsMultRPS"))->Fill(fFlowEvent->GetNumberOfRPs(),cc);
1552   }
1553   }
1554   // TPC vs VZERO
1555   ((TH3F*) fOutputFlowObs->FindObject( "hTPCVZE_AngleQ" ))->Fill(psi[0],psi[2],cc);
1556 }
1557
1558 //____________________________________________________________________________
1559 void AliAnalysisTaskSEHFQA::Terminate(Option_t */*option*/){
1560   //terminate analysis
1561
1562   fNEntries = dynamic_cast<TH1F*>(GetOutputData(1));
1563   if(!fNEntries){
1564     printf("ERROR: %s not available\n",GetOutputSlot(1)->GetContainer()->GetName());
1565     return;
1566   }
1567
1568   fOutputPID = dynamic_cast<TList*> (GetOutputData(2));
1569   if (!fOutputPID && fOnOff[1]) {
1570     printf("ERROR: %s not available\n",GetOutputSlot(2)->GetContainer()->GetName());
1571     return;
1572   }
1573
1574   fOutputTrack = dynamic_cast<TList*> (GetOutputData(3));
1575   if (!fOutputTrack && fOnOff[0]) {
1576     printf("ERROR: %s not available\n",GetOutputSlot(3)->GetContainer()->GetName());
1577     return;
1578   }
1579
1580 }
1581