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