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