PWG/FLOW/Tasks/AliAnalysisTaskVnV0.cxx

   1 #include "AliAnalysisTaskVnV0.h"
   2
   3 // ROOT includes
   4 #include <TMath.h>
   5
   6 // AliRoot includes
   7 #include "AliInputEventHandler.h"
   8 #include "AliAODEvent.h"
   9 #include "AliAODVertex.h"
  10 #include "AliAODTrack.h"
  11 #include "AliCentrality.h"
  12 #include "AliVHeader.h"
  13 #include "AliAODVZERO.h"
  14 #include "TFile.h"
  15 #include "AliOADBContainer.h"
  16 #include "TH2F.h"
  17 #include "TF1.h"
  18 #include "AliGenHijingEventHeader.h"
  19 #include "AliMCEvent.h"
  20 #include "AliAODMCHeader.h"
  21 #include "AliAODMCParticle.h"
  22 #include "TChain.h"
  23
  24 // STL includes
  25 //#include <iostream>
  26 //using namespace std;
  27
  28 ClassImp(AliAnalysisTaskVnV0)
  29
  30 //_____________________________________________________________________________
  31 AliAnalysisTaskVnV0::AliAnalysisTaskVnV0():
  32   AliAnalysisTaskSE(),
  33   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
  34   fEtaCut(0.8),   // cut on |eta| < fEtaCut
  35   fMinPt(0.15),   // cut on pt > fMinPt
  36   fV2(kTRUE),
  37   fV3(kTRUE),
  38   fIsMC(kFALSE),
  39   fQAsw(kFALSE),
  40   fRun(-1),
  41   fList(new TList()),
  42   fList2(new TList()),
  43   fList3(new TList()),
  44   fList4(new TList()),
  45   fMultV0(NULL),
  46   fV0Cpol(100),
  47   fV0Apol(100),
  48   fHResTPCv0A2(NULL),
  49   fHResTPCv0C2(NULL),
  50   fHResv0Cv0A2(NULL),
  51   fHResTPCv0A3(NULL),
  52   fHResTPCv0C3(NULL),
  53   fHResv0Cv0A3(NULL),
  54   fPhiRPv0A(NULL),
  55   fPhiRPv0C(NULL),
  56   fPhiRPv0Av3(NULL),
  57   fPhiRPv0Cv3(NULL),
  58   fQA(NULL),
  59   fQA2(NULL),
  60   fQAv3(NULL),
  61   fQA2v3(NULL),
  62   fPID(new AliFlowBayesianPID()),
  63   fTree(NULL),
  64   fCentrality(-1),
  65   evPlAngV0ACor2(0),
  66   evPlAngV0CCor2(0),
  67   evPlAng2(0),
  68   evPlAngV0ACor3(0),
  69   evPlAngV0CCor3(0),
  70   evPlAng3(0),
  71   fContAllChargesV0A(NULL),
  72   fContAllChargesV0C(NULL),
  73   fContAllChargesV0Av3(NULL),
  74   fContAllChargesV0Cv3(NULL),
  75   fContAllChargesMC(NULL),
  76   fHResMA2(NULL),
  77   fHResMC2(NULL),
  78   fHResAC2(NULL),
  79   fHResMA3(NULL),
  80   fHResMC3(NULL),
  81   fHResAC3(NULL),
  82   fContAllChargesMCA(NULL),
  83   fContAllChargesMCC(NULL),
  84   fContAllChargesMCAv3(NULL),
  85   fContAllChargesMCCv3(NULL),
  86   fFillDCA(kFALSE),
  87   fContQApid(NULL),
  88   fModulationDEDx(kFALSE)
  89 {
  90   // Default constructor (should not be used)
  91   fList->SetName("resultsV2");
  92   fList2->SetName("resultsV3");
  93   fList3->SetName("resultsMC");
  94   fList4->SetName("QA");
  95
  96   fList->SetOwner(kTRUE);
  97   fList2->SetOwner(kTRUE);
  98   fList3->SetOwner(kTRUE);
  99   fList4->SetOwner(kTRUE);
 100
 101   fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h
 102 }
 103
 104 //______________________________________________________________________________
 105 AliAnalysisTaskVnV0::AliAnalysisTaskVnV0(const char *name):
 106   AliAnalysisTaskSE(name),
 107   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
 108   fEtaCut(0.8),   // cut on |eta| < fEtaCut
 109   fMinPt(0.15),   // cut on pt > fMinPt
 110   fV2(kTRUE),
 111   fV3(kTRUE),
 112   fIsMC(kFALSE),
 113   fQAsw(kFALSE),
 114   fRun(-1),
 115   fList(new TList()),
 116   fList2(new TList()),
 117   fList3(new TList()),
 118   fList4(new TList()),
 119   fMultV0(NULL),
 120   fV0Cpol(100),
 121   fV0Apol(100),
 122   fHResTPCv0A2(NULL),
 123   fHResTPCv0C2(NULL),
 124   fHResv0Cv0A2(NULL),
 125   fHResTPCv0A3(NULL),
 126   fHResTPCv0C3(NULL),
 127   fHResv0Cv0A3(NULL),
 128   fPhiRPv0A(NULL),
 129   fPhiRPv0C(NULL),
 130   fPhiRPv0Av3(NULL),
 131   fPhiRPv0Cv3(NULL),
 132   fQA(NULL),
 133   fQA2(NULL),
 134   fQAv3(NULL),
 135   fQA2v3(NULL),
 136   fPID(new AliFlowBayesianPID()),
 137   fTree(NULL),
 138   fCentrality(-1),
 139   evPlAngV0ACor2(0),
 140   evPlAngV0CCor2(0),
 141   evPlAng2(0),
 142   evPlAngV0ACor3(0),
 143   evPlAngV0CCor3(0),
 144   evPlAng3(0),
 145   fContAllChargesV0A(NULL),
 146   fContAllChargesV0C(NULL),
 147   fContAllChargesV0Av3(NULL),
 148   fContAllChargesV0Cv3(NULL),
 149   fContAllChargesMC(NULL),
 150   fHResMA2(NULL),
 151   fHResMC2(NULL),
 152   fHResAC2(NULL),
 153   fHResMA3(NULL),
 154   fHResMC3(NULL),
 155   fHResAC3(NULL),
 156   fContAllChargesMCA(NULL),
 157   fContAllChargesMCC(NULL),
 158   fContAllChargesMCAv3(NULL),
 159   fContAllChargesMCCv3(NULL),
 160   fFillDCA(kFALSE),
 161   fContQApid(NULL),
 162   fModulationDEDx(kFALSE)
 163 {
 164
 165   DefineOutput(1, TList::Class());
 166   DefineOutput(2, TList::Class());
 167   DefineOutput(3, TList::Class());
 168   DefineOutput(4, TList::Class());
 169
 170   // Output slot #1 writes into a TTree
 171   fList->SetName("resultsV2");
 172   fList2->SetName("resultsV3");
 173   fList3->SetName("resultsMC");
 174   fList4->SetName("QA");
 175
 176   fList->SetOwner(kTRUE);
 177   fList2->SetOwner(kTRUE);
 178   fList3->SetOwner(kTRUE);
 179   fList4->SetOwner(kTRUE);
 180
 181   fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h
 182 }
 183 //_____________________________________________________________________________
 184 AliAnalysisTaskVnV0::~AliAnalysisTaskVnV0()
 185 {
 186
 187 }
 188
 189 //______________________________________________________________________________
 190 void AliAnalysisTaskVnV0::UserCreateOutputObjects()
 191 {
 192
 193   if(fIsMC) fPID->SetMC(kTRUE);
 194
 195
 196   // Tree for EP debug (comment the adding to v2 list id not needed)
 197   fTree = new TTree("tree","tree");
 198   fTree->Branch("evPlAngV0ACor2",&evPlAngV0ACor2,"evPlAngV0ACor2/F");
 199   fTree->Branch("evPlAngV0CCor2",&evPlAngV0CCor2,"evPlAngV0CCor2/F");
 200   fTree->Branch("evPlAng2",&evPlAng2,"evPlAng2/F");
 201   fTree->Branch("fCentrality",&fCentrality,"fCentrality/F");
 202   fTree->Branch("evPlAngV0ACor3",&evPlAngV0ACor3,"evPlAngV0ACor3/F");
 203   fTree->Branch("evPlAngV0CCor3",&evPlAngV0CCor3,"evPlAngV0CCor3/F");
 204   fTree->Branch("evPlAng3",&evPlAng3,"evPlAng3/F");
 205
 206
 207   // Container analyses (different steps mean different species)
 208   const Int_t nPtBinsTOF = 45;
 209   Double_t binsPtTOF[nPtBinsTOF+1] = {0., 0.05,  0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.25, 2.5, 2.75,3.0,3.25,3.5,3.75,4.0,4.5,5,5.5,6,6.5,7,8,9,10,12,15,20};
 210   const Int_t nCentrTOF = 9;
 211   const Int_t nPsiTOF = 10;
 212   const Int_t nChargeBinsTOFres = 2;
 213   const Int_t nCentrTOFres = 9;
 214   const Int_t nProbTOFres = 4;
 215   const Int_t nPsiTOFres = 10;
 216   const Int_t nMaskPID = 3;
 217
 218   Int_t nDCABin = 1; // put to 1 not to store this info
 219   if(fFillDCA) nDCABin = 3;
 220   if(fIsMC && nDCABin>1)  nDCABin = 6;
 221   /*
 222     0 = DCAxy < 2.4 && all (or Physical primary if MC)
 223     1 = DCAxy > 2.4 && all (or Physical primary if MC)
 224     2 = DCAxy < 2.4 && not Physical Primary for MC
 225     3 = DCAxy > 2.4 && not Physical Primary for MC
 226   */
 227
 228   Int_t binsTOF[6] = {nCentrTOFres,nChargeBinsTOFres,nProbTOFres,nPsiTOFres,nMaskPID,nDCABin};
 229   Int_t binsTOFmc[5] = {nCentrTOFres,nChargeBinsTOFres,1,nPsiTOFres,2};
 230   Int_t binsTOFmcPureMC[5] = {nCentrTOFres,nChargeBinsTOFres,1,nPsiTOFres,1};
 231
 232   // v2 container
 233   fContAllChargesV0A = new AliFlowVZEROResults("v2A",6,binsTOF);
 234   fContAllChargesV0A->SetVarRange(0,-0.5,8.5); // centrality
 235   fContAllChargesV0A->SetVarRange(1,-1.5,1.5);  // charge
 236   fContAllChargesV0A->SetVarRange(2,0.6,1.0001);// prob
 237   fContAllChargesV0A->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
 238   fContAllChargesV0A->SetVarRange(4,-0.5,2.5); // pid mask
 239   fContAllChargesV0A->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
 240   fContAllChargesV0A->SetVarName(0,"centrality");
 241   fContAllChargesV0A->SetVarName(1,"charge");
 242   fContAllChargesV0A->SetVarName(2,"prob");
 243   fContAllChargesV0A->SetVarName(3,"#Psi");
 244   fContAllChargesV0A->SetVarName(4,"PIDmask");
 245   fContAllChargesV0A->SetVarName(5,"DCAbin");
 246   if(fV2) fContAllChargesV0A->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 247   if(fV2) fContAllChargesV0A->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 248   if(fV2) fContAllChargesV0A->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 249   if(fV2) fContAllChargesV0A->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 250   if(fV2) fContAllChargesV0A->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 251   if(fV2) fContAllChargesV0A->AddSpecies("d",nPtBinsTOF,binsPtTOF);
 252   if(fV2) fContAllChargesV0A->AddSpecies("t",nPtBinsTOF,binsPtTOF);
 253   if(fV2) fContAllChargesV0A->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
 254   if(fV2) fContAllChargesV0A->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 255
 256   fContAllChargesV0C = new AliFlowVZEROResults("v2C",6,binsTOF);
 257   fContAllChargesV0C->SetVarRange(0,-0.5,8.5); // centrality
 258   fContAllChargesV0C->SetVarRange(1,-1.5,1.5);  // charge
 259   fContAllChargesV0C->SetVarRange(2,0.6,1.0001);// prob
 260   fContAllChargesV0C->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
 261   fContAllChargesV0C->SetVarRange(4,-0.5,2.5); // pid mask
 262   fContAllChargesV0C->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
 263   fContAllChargesV0C->SetVarName(0,"centrality");
 264   fContAllChargesV0C->SetVarName(1,"charge");
 265   fContAllChargesV0C->SetVarName(2,"prob");
 266   fContAllChargesV0C->SetVarName(3,"#Psi");
 267   fContAllChargesV0C->SetVarName(4,"PIDmask");
 268   fContAllChargesV0C->SetVarName(5,"DCAbin");
 269   if(fV2) fContAllChargesV0C->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 270   if(fV2) fContAllChargesV0C->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 271   if(fV2) fContAllChargesV0C->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 272   if(fV2) fContAllChargesV0C->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 273   if(fV2) fContAllChargesV0C->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 274   if(fV2) fContAllChargesV0C->AddSpecies("d",nPtBinsTOF,binsPtTOF);
 275   if(fV2) fContAllChargesV0C->AddSpecies("t",nPtBinsTOF,binsPtTOF);
 276   if(fV2) fContAllChargesV0C->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
 277   if(fV2) fContAllChargesV0C->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 278
 279   fList->Add(fContAllChargesV0A);
 280   fList->Add(fContAllChargesV0C);
 281
 282   if(fIsMC && fV2){
 283     fContAllChargesMC = new AliFlowVZEROResults("v2mc",5,binsTOFmc);
 284     fContAllChargesMC->SetVarRange(0,-0.5,8.5); // centrality
 285     fContAllChargesMC->SetVarRange(1,-1.5,1.5);  // charge
 286     fContAllChargesMC->SetVarRange(2,0.6,1.0001);// prob
 287     fContAllChargesMC->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
 288     fContAllChargesMC->SetVarRange(4,-0.5,1.5); // pid mask
 289     fContAllChargesMC->SetVarName(0,"centrality");
 290     fContAllChargesMC->SetVarName(1,"charge");
 291     fContAllChargesMC->SetVarName(2,"prob");
 292     fContAllChargesMC->SetVarName(3,"#Psi");
 293     fContAllChargesMC->SetVarName(4,"PIDmask");
 294     fContAllChargesMC->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 295     fContAllChargesMC->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 296     fContAllChargesMC->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 297     fContAllChargesMC->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 298     fContAllChargesMC->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 299     fContAllChargesMC->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 300     fList3->Add(fContAllChargesMC);
 301
 302     fContAllChargesMCA = new AliFlowVZEROResults("v2mcA",5,binsTOFmcPureMC);
 303     fContAllChargesMCA->SetVarRange(0,-0.5,8.5); // centrality
 304     fContAllChargesMCA->SetVarRange(1,-1.5,1.5);  // charge
 305     fContAllChargesMCA->SetVarRange(2,0.6,1.0001);// prob
 306     fContAllChargesMCA->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
 307     fContAllChargesMCA->SetVarRange(4,-0.5,1.5); // pid mask
 308     fContAllChargesMCA->SetVarName(0,"centrality");
 309     fContAllChargesMCA->SetVarName(1,"charge");
 310     fContAllChargesMCA->SetVarName(2,"prob");
 311     fContAllChargesMCA->SetVarName(3,"#Psi");
 312     fContAllChargesMCA->SetVarName(4,"PIDmask");
 313     fContAllChargesMCA->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 314     fContAllChargesMCA->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 315     fContAllChargesMCA->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 316     fContAllChargesMCA->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 317     fContAllChargesMCA->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 318     fContAllChargesMCA->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 319     fList3->Add(fContAllChargesMCA);
 320
 321     fContAllChargesMCC = new AliFlowVZEROResults("v2mcC",5,binsTOFmcPureMC);
 322     fContAllChargesMCC->SetVarRange(0,-0.5,8.5); // centrality
 323     fContAllChargesMCC->SetVarRange(1,-1.5,1.5);  // charge
 324     fContAllChargesMCC->SetVarRange(2,0.6,1.0001);// prob
 325     fContAllChargesMCC->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
 326     fContAllChargesMCC->SetVarRange(4,-0.5,1.5); // pid mask
 327     fContAllChargesMCC->SetVarName(0,"centrality");
 328     fContAllChargesMCC->SetVarName(1,"charge");
 329     fContAllChargesMCC->SetVarName(2,"prob");
 330     fContAllChargesMCC->SetVarName(3,"#Psi");
 331     fContAllChargesMCC->SetVarName(4,"PIDmask");
 332     fContAllChargesMCC->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 333     fContAllChargesMCC->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 334     fContAllChargesMCC->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 335     fContAllChargesMCC->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 336     fContAllChargesMCC->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 337     fContAllChargesMCC->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 338     fList3->Add(fContAllChargesMCC);
 339   }
 340
 341   // v3 container
 342   fContAllChargesV0Av3 = new AliFlowVZEROResults("v3A",6,binsTOF);
 343   fContAllChargesV0Av3->SetVarRange(0,-0.5,8.5); // centrality
 344   fContAllChargesV0Av3->SetVarRange(1,-1.5,1.5);  // charge
 345   fContAllChargesV0Av3->SetVarRange(2,0.6,1.0001);// prob
 346   fContAllChargesV0Av3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
 347   fContAllChargesV0Av3->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
 348   fContAllChargesV0Av3->SetVarRange(4,-0.5,2.5); // pid mask
 349   fContAllChargesV0Av3->SetVarName(0,"centrality");
 350   fContAllChargesV0Av3->SetVarName(1,"charge");
 351   fContAllChargesV0Av3->SetVarName(2,"prob");
 352   fContAllChargesV0Av3->SetVarName(3,"#Psi");
 353   fContAllChargesV0Av3->SetVarName(4,"PIDmask");
 354   fContAllChargesV0Av3->SetVarName(5,"DCAbin");
 355   if(fV3) fContAllChargesV0Av3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 356   if(fV3) fContAllChargesV0Av3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 357   if(fV3) fContAllChargesV0Av3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 358   if(fV3) fContAllChargesV0Av3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 359   if(fV3) fContAllChargesV0Av3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 360   if(fV3) fContAllChargesV0Av3->AddSpecies("d",nPtBinsTOF,binsPtTOF);
 361   if(fV3) fContAllChargesV0Av3->AddSpecies("t",nPtBinsTOF,binsPtTOF);
 362   if(fV3) fContAllChargesV0Av3->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
 363   if(fV3) fContAllChargesV0Av3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 364
 365   fContAllChargesV0Cv3 = new AliFlowVZEROResults("v3C",6,binsTOF);
 366   fContAllChargesV0Cv3->SetVarRange(0,-0.5,8.5); // centrality
 367   fContAllChargesV0Cv3->SetVarRange(1,-1.5,1.5);  // charge
 368   fContAllChargesV0Cv3->SetVarRange(2,0.6,1.0001);// prob
 369   fContAllChargesV0Cv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
 370   fContAllChargesV0Cv3->SetVarRange(4,-0.5,2.5); // pid mask
 371   fContAllChargesV0Cv3->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
 372   fContAllChargesV0Cv3->SetVarName(0,"centrality");
 373   fContAllChargesV0Cv3->SetVarName(1,"charge");
 374   fContAllChargesV0Cv3->SetVarName(2,"prob");
 375   fContAllChargesV0Cv3->SetVarName(3,"#Psi");
 376   fContAllChargesV0Cv3->SetVarName(4,"PIDmask");
 377   fContAllChargesV0Cv3->SetVarName(5,"DCAbin");
 378   if(fV3) fContAllChargesV0Cv3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 379   if(fV3) fContAllChargesV0Cv3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 380   if(fV3) fContAllChargesV0Cv3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 381   if(fV3) fContAllChargesV0Cv3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 382   if(fV3) fContAllChargesV0Cv3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 383   if(fV3) fContAllChargesV0Cv3->AddSpecies("d",nPtBinsTOF,binsPtTOF);
 384   if(fV3) fContAllChargesV0Cv3->AddSpecies("t",nPtBinsTOF,binsPtTOF);
 385   if(fV3) fContAllChargesV0Cv3->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
 386   if(fV3) fContAllChargesV0Cv3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 387
 388   fList2->Add(fContAllChargesV0Av3);
 389   fList2->Add(fContAllChargesV0Cv3);
 390
 391   if(fIsMC && fV3){
 392     fContAllChargesMCAv3 = new AliFlowVZEROResults("v3mcA",5,binsTOFmcPureMC);
 393     fContAllChargesMCAv3->SetVarRange(0,-0.5,8.5); // centrality
 394     fContAllChargesMCAv3->SetVarRange(1,-1.5,1.5);  // charge
 395     fContAllChargesMCAv3->SetVarRange(2,0.6,1.0001);// prob
 396     fContAllChargesMCAv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
 397     fContAllChargesMCAv3->SetVarRange(4,-0.5,1.5); // pid mask
 398     fContAllChargesMCAv3->SetVarName(0,"centrality");
 399     fContAllChargesMCAv3->SetVarName(1,"charge");
 400     fContAllChargesMCAv3->SetVarName(2,"prob");
 401     fContAllChargesMCAv3->SetVarName(3,"#Psi");
 402     fContAllChargesMCAv3->SetVarName(4,"PIDmask");
 403     fContAllChargesMCAv3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 404     fContAllChargesMCAv3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 405     fContAllChargesMCAv3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 406     fContAllChargesMCAv3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 407     fContAllChargesMCAv3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 408     fContAllChargesMCAv3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 409     fList3->Add(fContAllChargesMCAv3);
 410
 411     fContAllChargesMCCv3 = new AliFlowVZEROResults("v3mcC",5,binsTOFmcPureMC);
 412     fContAllChargesMCCv3->SetVarRange(0,-0.5,8.5); // centrality
 413     fContAllChargesMCCv3->SetVarRange(1,-1.5,1.5);  // charge
 414     fContAllChargesMCCv3->SetVarRange(2,0.6,1.0001);// prob
 415     fContAllChargesMCCv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
 416     fContAllChargesMCCv3->SetVarRange(4,-0.5,1.5); // pid mask
 417     fContAllChargesMCCv3->SetVarName(0,"centrality");
 418     fContAllChargesMCCv3->SetVarName(1,"charge");
 419     fContAllChargesMCCv3->SetVarName(2,"prob");
 420     fContAllChargesMCCv3->SetVarName(3,"#Psi");
 421     fContAllChargesMCCv3->SetVarName(4,"PIDmask");
 422     fContAllChargesMCCv3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 423     fContAllChargesMCCv3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 424     fContAllChargesMCCv3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 425     fContAllChargesMCCv3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 426     fContAllChargesMCCv3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 427     fContAllChargesMCCv3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 428     fList3->Add(fContAllChargesMCCv3);
 429   }
 430
 431   // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
 432   // v2
 433   fHResTPCv0A2 = new TProfile("hResTPCv0A2","",9,0,9);
 434   fHResTPCv0C2 = new TProfile("hResTPCv0C2","",9,0,9);
 435   fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",9,0,9);
 436
 437   fList->Add(fHResTPCv0A2);
 438   fList->Add(fHResTPCv0C2);
 439   fList->Add(fHResv0Cv0A2);
 440
 441   // v3
 442   fHResTPCv0A3 = new TProfile("hResTPCv0A3","",9,0,9);
 443   fHResTPCv0C3 = new TProfile("hResTPCv0C3","",9,0,9);
 444   fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",9,0,9);
 445
 446   fList2->Add(fHResTPCv0A3);
 447   fList2->Add(fHResTPCv0C3);
 448   fList2->Add(fHResv0Cv0A3);
 449
 450   // MC as in the dataEP resolution (but using MC tracks)
 451   if(fIsMC && fV3){
 452     fHResMA2 = new TProfile("hResMA2","",9,0,9);
 453     fHResMC2 = new TProfile("hResMC2","",9,0,9);
 454     fHResAC2 = new TProfile("hResAC2","",9,0,9);
 455     fList3->Add(fHResMA2);
 456     fList3->Add(fHResMC2);
 457     fList3->Add(fHResAC2);
 458   }
 459   if(fIsMC && fV3){
 460     fHResMA3 = new TProfile("hResMA3","",9,0,9);
 461     fHResMC3 = new TProfile("hResMC3","",9,0,9);
 462     fHResAC3 = new TProfile("hResAC3","",9,0,9);
 463     fList3->Add(fHResMA3);
 464     fList3->Add(fHResMC3);
 465     fList3->Add(fHResAC3);
 466   }
 467
 468
 469   // V0A and V0C event plane distributions
 470   //v2
 471   fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
 472   fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
 473
 474   //v3
 475   fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
 476   fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
 477
 478   // QA container
 479   // v2
 480   const Int_t nDETsignal = 50;
 481   Double_t binDETsignal[nDETsignal+1];
 482   for(Int_t i=0;i<nDETsignal+1;i++){
 483     binDETsignal[i] = -5 + i*10. / nDETsignal;
 484   }
 485 //   const Int_t nEta = 5;
 486 //   Double_t binEta[nEta+1];
 487 //   for(Int_t i=0;i<nEta+1;i++){
 488 //     binEta[i] = -1 + i*2. / nEta;
 489 //   }
 490
 491   const Int_t nDeltaPhi = 5;
 492   const Int_t nDeltaPhiV3 = 7;
 493
 494   Int_t binsQA[5] = {nCentrTOF,7,5,nDeltaPhi,2};
 495   Int_t binsQAv3[5] = {nCentrTOF,7,5,nDeltaPhiV3,2};
 496
 497
 498   fQA = new AliFlowVZEROQA("v2AQA",5,binsQA);
 499   fQA->SetVarRange(0,-0.5,8.5); // centrality
 500   fQA->SetVarRange(1,0,7);  // pt
 501   fQA->SetVarRange(2,0.,1.0001);// prob
 502   fQA->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
 503   fQA->SetVarRange(4,-0.5,1.5); // pid mask
 504   fQA->SetVarName(0,"centrality");
 505   fQA->SetVarName(1,"p_{t}");
 506   fQA->SetVarName(2,"prob");
 507   fQA->SetVarName(3,"#Psi");
 508   fQA->SetVarName(4,"PIDmask");
 509   if(fQAsw && fV2) fQA->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 510   if(fQAsw && fV2) fQA->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 511   if(fQAsw && fV2) fQA->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 512 //   if(fQAsw && fV2) fQA->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 513 //   fQA->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 514 //   fQA->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 515 //   fQA->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 516
 517   fQA2 = new AliFlowVZEROQA("v2CQA",5,binsQA);
 518   fQA2->SetVarRange(0,-0.5,8.5); // centrality
 519   fQA2->SetVarRange(1,0,7);  // pt
 520   fQA2->SetVarRange(2,0.,1.0001);// prob
 521   fQA2->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
 522   fQA2->SetVarRange(4,-0.5,1.5); // pid mask
 523   fQA2->SetVarName(0,"centrality");
 524   fQA2->SetVarName(1,"p_{t}");
 525   fQA2->SetVarName(2,"prob");
 526   fQA2->SetVarName(3,"#Psi");
 527   fQA2->SetVarName(4,"PIDmask");
 528   if(fQAsw && fV2) fQA2->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 529   if(fQAsw && fV2) fQA2->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 530   if(fQAsw && fV2) fQA2->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 531 //   if(fQAsw && fV2) fQA2->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 532 //   fQA2->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 533 //   fQA2->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 534 //   fQA2->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 535
 536   fQAv3 = new AliFlowVZEROQA("v3AQA",5,binsQAv3);
 537   fQAv3->SetVarRange(0,-0.5,8.5); // centrality
 538   fQAv3->SetVarRange(1,0,7);  // pt
 539   fQAv3->SetVarRange(2,0.,1.0001);// prob
 540   fQAv3->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
 541   fQAv3->SetVarRange(4,-0.5,1.5); // pid mask
 542   fQAv3->SetVarName(0,"centrality");
 543   fQAv3->SetVarName(1,"p_{t}");
 544   fQAv3->SetVarName(2,"prob");
 545   fQAv3->SetVarName(3,"#Psi");
 546   fQAv3->SetVarName(4,"PIDmask");
 547   if(fQAsw && fV3) fQAv3->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 548 //   if(fQAsw && fV3) fQAv3->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 549 //   if(fQAsw && fV3) fQAv3->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 550 //   if(fQAsw && fV2) fQAv3->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 551 //   fQAv3->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 552 //   fQAv3->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 553 //   fQAv3->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 554
 555   fQA2v3 = new AliFlowVZEROQA("v3CQA",5,binsQAv3);
 556   fQA2v3->SetVarRange(0,-0.5,8.5); // centrality
 557   fQA2v3->SetVarRange(1,0,7);  // pt
 558   fQA2v3->SetVarRange(2,0.,1.0001);// prob
 559   fQA2v3->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
 560   fQA2v3->SetVarRange(4,-0.5,1.5); // pid mask
 561   fQA2v3->SetVarName(0,"centrality");
 562   fQA2v3->SetVarName(1,"p_{t}");
 563   fQA2v3->SetVarName(2,"prob");
 564   fQA2v3->SetVarName(3,"#Psi");
 565   fQA2v3->SetVarName(4,"PIDmask");
 566   if(fQAsw && fV3) fQA2v3->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 567 //   if(fQAsw && fV3) fQA2v3->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 568 //   if(fQAsw && fV3) fQA2v3->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 569 //   if(fQAsw && fV2) fQA2v3->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 570 //   fQA2v3->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 571 //   fQA2v3->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 572 //   fQA2v3->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 573
 574   fList->Add(fPhiRPv0A);
 575   fList->Add(fPhiRPv0C);
 576
 577   if(fQAsw && fV2){
 578     fList4->Add(fQA);
 579     fList4->Add(fQA2);
 580   }
 581
 582   fList2->Add(fPhiRPv0Av3);
 583   fList2->Add(fPhiRPv0Cv3);
 584
 585   if(fQAsw && fV3){
 586    fList4->Add(fQAv3);
 587    fList4->Add(fQA2v3);
 588   }
 589
 590   //  fList->Add(fTree); // comment if not needed
 591
 592   const Int_t nDCA = 300;
 593   Double_t DCAbin[nDCA+1];
 594   for(Int_t i=0;i <= nDCA;i++){
 595     DCAbin[i] = -3 +i*6.0/nDCA;
 596   }
 597
 598   char nameHistos[100];
 599   for(Int_t iC=0;iC < nCentrBin;iC++){
 600     snprintf(nameHistos,100,"fHdcaPtPiCent%i",iC);
 601     fHdcaPt[iC][0] = new TH2D(nameHistos,"DCA_{xy} for #pi;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 602     snprintf(nameHistos,100,"fHdcaPtKaCent%i",iC);
 603     fHdcaPt[iC][1] = new TH2D(nameHistos,"DCA_{xy} for K;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 604     snprintf(nameHistos,100,"fHdcaPtPrCent%i",iC);
 605     fHdcaPt[iC][2] = new TH2D(nameHistos,"DCA_{xy} for #bar{p};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 606     snprintf(nameHistos,100,"fHdcaPtElCent%i",iC);
 607     fHdcaPt[iC][3] = new TH2D(nameHistos,"DCA_{xy} for e;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 608     snprintf(nameHistos,100,"fHdcaPtDeCent%i",iC);
 609     fHdcaPt[iC][4] = new TH2D(nameHistos,"DCA_{xy} for #bar{d};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 610     snprintf(nameHistos,100,"fHdcaPtTrCent%i",iC);
 611     fHdcaPt[iC][5] = new TH2D(nameHistos,"DCA_{xy} for #bar{t};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 612     snprintf(nameHistos,100,"fHdcaPtHeCent%i",iC);
 613     fHdcaPt[iC][6] = new TH2D(nameHistos,"DCA_{xy} for #bar{^{3}He};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 614   }
 615
 616   if(fFillDCA && fQAsw){
 617     for(Int_t i=0;i<7;i++)
 618       for(Int_t iC=0;iC < nCentrBin;iC++)
 619         fList4->Add(fHdcaPt[iC][i]);
 620   }
 621   if(fIsMC){
 622     for(Int_t iC=0;iC < nCentrBin;iC++){
 623       snprintf(nameHistos,100,"fHdcaPtPiSecCent%i",iC);
 624       fHdcaPtSec[iC][0] = new TH2D(nameHistos,"DCA_{xy} for secondary #pi;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 625       snprintf(nameHistos,100,"fHdcaPtKaSecCent%i",iC);
 626       fHdcaPtSec[iC][1] = new TH2D(nameHistos,"DCA_{xy} for secondary K;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 627       snprintf(nameHistos,100,"fHdcaPtPrSecCent%i",iC);
 628       fHdcaPtSec[iC][2] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{p};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 629       snprintf(nameHistos,100,"fHdcaPtElSecCent%i",iC);
 630       fHdcaPtSec[iC][3] = new TH2D(nameHistos,"DCA_{xy} for secondary e;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 631       snprintf(nameHistos,100,"fHdcaPtDeSecCent%i",iC);
 632       fHdcaPtSec[iC][4] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{d};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 633       snprintf(nameHistos,100,"fHdcaPtTrSecCent%i",iC);
 634       fHdcaPtSec[iC][5] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{t};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 635       snprintf(nameHistos,100,"fHdcaPtHeSecCent%i",iC);
 636       fHdcaPtSec[iC][6] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{^{3}He};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 637     }
 638
 639     if(fFillDCA && fQAsw){
 640       for(Int_t i=0;i<7;i++)
 641         for(Int_t iC=0;iC < nCentrBin;iC++)
 642           fList4->Add(fHdcaPtSec[iC][i]);
 643     }
 644   }
 645
 646   // Add TProfile Extra QA
 647   const Int_t nBinQApid = 2;
 648   Int_t binQApid[nBinQApid] = {nCentrTOF,200};
 649   const Int_t nbinsigma = 100;
 650   Double_t nsigmaQA[nbinsigma+1];
 651   for(Int_t i=0;i<nbinsigma+1;i++){
 652     nsigmaQA[i] = -10 + 20.0*i/nbinsigma;
 653   }
 654   fContQApid = new AliFlowVZEROResults("qaPID",nBinQApid,binQApid);
 655   fContQApid->SetVarRange(0,-0.5,8.5); // centrality
 656   fContQApid->SetVarRange(1,0,20);  // charge
 657   fContQApid->SetVarName(0,"centrality");
 658   fContQApid->SetVarName(1,"p_{t}");
 659   fContQApid->AddSpecies("piTPC",nbinsigma,nsigmaQA);
 660   fContQApid->AddSpecies("piTOF",nbinsigma,nsigmaQA);
 661   fContQApid->AddSpecies("kaTPC",nbinsigma,nsigmaQA);
 662   fContQApid->AddSpecies("kaTOF",nbinsigma,nsigmaQA);
 663   fContQApid->AddSpecies("prTPC",nbinsigma,nsigmaQA);
 664   fContQApid->AddSpecies("prTOF",nbinsigma,nsigmaQA);
 665   if(fV2) fList->Add(fContQApid);
 666
 667   printf("Output creation ok!!\n\n\n\n");
 668
 669   // Post output data.
 670   if(fV2) PostData(1, fList);
 671   if(fV3) PostData(2, fList2);
 672   if(fIsMC) PostData(3, fList3);
 673   if(fQAsw) PostData(4, fList4);
 674
 675 }
 676
 677 //______________________________________________________________________________
 678 void AliAnalysisTaskVnV0::UserExec(Option_t *)
 679 {
 680     // Main loop
 681     // Called for each event
 682
 683     fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 684     if(!fOutputAOD){
 685         Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 686         this->Dump();
 687         return;
 688     }
 689
 690     Int_t run = fOutputAOD->GetRunNumber();
 691
 692     if(run != fRun){
 693         // Load the calibrations run dependent
 694         OpenInfoCalbration(run);
 695         fRun=run;
 696     }
 697
 698     Float_t zvtx = GetVertex(fOutputAOD);
 699
 700
 701
 702     //Get the MC object
 703     if(fIsMC){
 704       AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 705       if (!mcHeader) {
 706         AliError("Could not find MC Header in AOD");
 707         return;
 708       }
 709     }
 710
 711     /*
 712     AliMCEvent* mcEvent = MCEvent();
 713     if (!mcEvent) {
 714       Printf("ERROR: Could not retrieve MC event");
 715       return;
 716     }
 717
 718     Double_t gReactionPlane = -999., gImpactParameter = -999.;
 719     //Get the MC header
 720     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
 721     if (headerH) {
 722       //Printf("=====================================================");
 723       //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());
 724       //Printf("=====================================================");
 725       gReactionPlane = headerH->ReactionPlaneAngle();
 726       gImpactParameter = headerH->ImpactParameter();
 727     }
 728
 729 */
 730
 731     if (TMath::Abs(zvtx) < fVtxCut) {
 732       //Centrality
 733       Float_t v0Centr  = -10.;
 734       Float_t trkCentr  = -10.;
 735       AliCentrality *centrality = fOutputAOD->GetCentrality();
 736       if (centrality){
 737         v0Centr  = centrality->GetCentralityPercentile("V0M");
 738         trkCentr = centrality->GetCentralityPercentile("TRK");
 739       }
 740
 741       if(TMath::Abs(v0Centr - trkCentr) < 5.0){ // consistency cut on centrality selection
 742         fPID->SetDetResponse(fOutputAOD, v0Centr); // Set the PID object for each event!!!!
 743         Analyze(fOutputAOD,v0Centr); // Do analysis!!!!
 744
 745         fCentrality = v0Centr;
 746         if(fV2) fTree->Fill();
 747       }
 748     }
 749
 750 }
 751
 752 //________________________________________________________________________
 753 void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
 754 {
 755   Float_t mass[8] = {5.10998909999999971e-04, 1.05658000000000002e-01, 1.39570000000000000e-01, 4.93676999999999977e-01, 9.38271999999999995e-01,1.87783699999999998,2.81740199999999996,1.40805449999999999};
 756
 757   // Event plane resolution for v2
 758   Float_t evPlRes[18] = {0.350582,0.505393,0.607845,0.632913,0.592230,0.502489,0.381717,0.249539,0.133180, // V0A vs. centrality
 759                          0.446480,0.612705,0.712222,0.736200,0.697907,0.610114,0.481009,0.327402,0.182277};// V0C vs. centrality
 760
 761   Int_t iC = -1;
 762   if (v0Centr >0 && v0Centr < 80){ // analysis only for 0-80% centrality classes
 763     // centrality bins
 764     if(v0Centr < 5) iC = 0;
 765     else if(v0Centr < 10) iC = 1;
 766     else if(v0Centr < 20) iC = 2;
 767     else if(v0Centr < 30) iC = 3;
 768     else if(v0Centr < 40) iC = 4;
 769     else if(v0Centr < 50) iC = 5;
 770     else if(v0Centr < 60) iC = 6;
 771     else if(v0Centr < 70) iC = 7;
 772     else iC = 8;
 773
 774     //reset Q vector info
 775     Double_t Qxa2 = 0, Qya2 = 0;
 776     Double_t Qxc2 = 0, Qyc2 = 0;
 777     Double_t Qxa3 = 0, Qya3 = 0;
 778     Double_t Qxc3 = 0, Qyc3 = 0;
 779
 780     Int_t nAODTracks = aodEvent->GetNumberOfTracks();
 781
 782     AliAODMCHeader *mcHeader = NULL;
 783     TClonesArray *mcArray = NULL;
 784     Float_t evplaneMC = 0;
 785     if(fIsMC){
 786       mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 787
 788       if (mcHeader) {
 789         evplaneMC = mcHeader->GetReactionPlaneAngle();
 790         if(evplaneMC > TMath::Pi()/2 && evplaneMC <=  TMath::Pi()*3/2) evplaneMC-=TMath::Pi();
 791         else if(evplaneMC > TMath::Pi()*3/2) evplaneMC-=2*TMath::Pi();
 792         mcArray = (TClonesArray*)fOutputAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 793
 794         if(mcArray){
 795           Float_t QxMCv2[3] = {0,0,0};
 796           Float_t QyMCv2[3] = {0,0,0};
 797           Float_t QxMCv3[3] = {0,0,0};
 798           Float_t QyMCv3[3] = {0,0,0};
 799           Float_t EvPlaneMCV2[3] = {0,0,0};
 800           Float_t EvPlaneMCV3[3] = {0,0,0};
 801           Float_t etaMin[3] = {2.8,-3.6,-0.8}; // A-side, C-side M-barrel
 802           Float_t etaMax[3] = {4.88,-1.8,0.8};
 803
 804           // analysis on MC tracks
 805           Int_t nMCtrack = mcArray->GetEntries() ;
 806
 807           // EP computation with MC tracks
 808           for(Int_t iT=0;iT < nMCtrack;iT++){
 809             AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
 810             if(!mctr || !(mctr->IsPrimary()) || !(mctr->Charge()) || mctr->Pt() < 0.2) continue;
 811
 812             Float_t eta = mctr->Eta();
 813
 814             for(Int_t iD=0;iD<3;iD++){
 815               if(eta > etaMin[iD] && eta < etaMax[iD]){
 816                 Float_t phi = mctr->Phi();
 817                 QxMCv2[iD] += TMath::Cos(2*phi);
 818                 QyMCv2[iD] += TMath::Sin(2*phi);
 819                 QxMCv3[iD] += TMath::Cos(3*phi);
 820                 QyMCv3[iD] += TMath::Sin(3*phi);
 821               }
 822             }
 823           }
 824
 825           if(fV2){
 826             EvPlaneMCV2[0] = TMath::ATan2(QyMCv2[0],QxMCv2[0])/2.;
 827             EvPlaneMCV2[1] = TMath::ATan2(QyMCv2[1],QxMCv2[1])/2.;
 828             EvPlaneMCV2[2] = TMath::ATan2(QyMCv2[2],QxMCv2[2])/2.;
 829             fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
 830             fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
 831             fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
 832           }
 833           if(fV3){
 834             EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
 835             EvPlaneMCV3[1] = TMath::ATan2(QyMCv3[1],QxMCv3[1])/3.;
 836             EvPlaneMCV3[2] = TMath::ATan2(QyMCv3[2],QxMCv3[2])/3.;
 837             fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
 838             fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
 839             fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
 840           }
 841
 842           // flow A and C side
 843           Float_t xMCepAv2[5] = {iC,0/*charge*/,1,EvPlaneMCV2[0],1};
 844           Float_t xMCepCv2[5] = {iC,0/*charge*/,1,EvPlaneMCV2[1],1};
 845           Float_t xMCepAv3[5] = {iC,0/*charge*/,1,EvPlaneMCV3[0],1};
 846           Float_t xMCepCv3[5] = {iC,0/*charge*/,1,EvPlaneMCV3[1],1};
 847
 848           for(Int_t iT=0;iT < nMCtrack;iT++){
 849             AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
 850             if(!mctr || !(mctr->IsPhysicalPrimary()) || !(mctr->Charge()) || TMath::Abs(mctr->Eta()) > 0.8 || mctr->Pt() < 0.2) continue;
 851             Int_t iS = TMath::Abs(mctr->GetPdgCode());
 852             Int_t charge = mctr->Charge();
 853             Float_t pt = mctr->Pt();
 854             Float_t phi = mctr->Phi();
 855
 856             if(charge > 0){
 857               xMCepAv2[1] = 1;
 858               xMCepCv2[1] = 1;
 859               xMCepAv3[1] = 1;
 860               xMCepCv3[1] = 1;
 861             }
 862             else{
 863               xMCepAv2[1] = -1;
 864               xMCepCv2[1] = -1;
 865               xMCepAv3[1] = -1;
 866               xMCepCv3[1] = -1;
 867             }
 868
 869             fContAllChargesMCA->Fill(0,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
 870             fContAllChargesMCC->Fill(0,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
 871             fContAllChargesMCAv3->Fill(0,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
 872             fContAllChargesMCCv3->Fill(0,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
 873
 874             if(iS==11){
 875               fContAllChargesMCA->Fill(4,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
 876               fContAllChargesMCC->Fill(4,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
 877               fContAllChargesMCAv3->Fill(4,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
 878               fContAllChargesMCCv3->Fill(4,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
 879             }
 880             else if(iS==13){
 881               fContAllChargesMCA->Fill(5,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
 882               fContAllChargesMCC->Fill(5,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
 883               fContAllChargesMCAv3->Fill(5,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
 884               fContAllChargesMCCv3->Fill(5,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
 885             }
 886             else if(iS==211){
 887               fContAllChargesMCA->Fill(1,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
 888               fContAllChargesMCC->Fill(1,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
 889               fContAllChargesMCAv3->Fill(1,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
 890               fContAllChargesMCCv3->Fill(1,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
 891             }
 892             else if(iS==321){
 893               fContAllChargesMCA->Fill(2,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
 894               fContAllChargesMCC->Fill(2,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
 895               fContAllChargesMCAv3->Fill(2,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
 896               fContAllChargesMCCv3->Fill(2,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
 897             }
 898             else if(iS==2212){
 899               fContAllChargesMCA->Fill(3,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
 900               fContAllChargesMCC->Fill(3,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
 901               fContAllChargesMCAv3->Fill(3,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
 902               fContAllChargesMCCv3->Fill(3,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
 903             }
 904           }
 905         }
 906       }
 907     }
 908
 909     //V0 info
 910     AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
 911
 912     for (Int_t iv0 = 0; iv0 < 64; iv0++) {
 913       Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
 914       Float_t multv0 = aodV0->GetMultiplicity(iv0);
 915       if (iv0 < 32){ // V0C
 916         Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
 917         Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
 918         Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
 919         Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
 920       } else {       // V0A
 921         Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
 922         Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
 923         Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
 924         Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
 925       }
 926     }
 927
 928     //grab for each centrality the proper histo with the Qx and Qy to do the recentering
 929     Double_t Qxamean2 = fMeanQ[iC][1][0];
 930     Double_t Qxarms2  = fWidthQ[iC][1][0];
 931     Double_t Qyamean2 = fMeanQ[iC][1][1];
 932     Double_t Qyarms2  = fWidthQ[iC][1][1];
 933     Double_t Qxamean3 = fMeanQv3[iC][1][0];
 934     Double_t Qxarms3  = fWidthQv3[iC][1][0];
 935     Double_t Qyamean3 = fMeanQv3[iC][1][1];
 936     Double_t Qyarms3  = fWidthQv3[iC][1][1];
 937
 938     Double_t Qxcmean2 = fMeanQ[iC][0][0];
 939     Double_t Qxcrms2  = fWidthQ[iC][0][0];
 940     Double_t Qycmean2 = fMeanQ[iC][0][1];
 941     Double_t Qycrms2  = fWidthQ[iC][0][1];
 942     Double_t Qxcmean3 = fMeanQv3[iC][0][0];
 943     Double_t Qxcrms3  = fWidthQv3[iC][0][0];
 944     Double_t Qycmean3 = fMeanQv3[iC][0][1];
 945     Double_t Qycrms3  = fWidthQv3[iC][0][1];
 946
 947     Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
 948     Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
 949     Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
 950     Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
 951     Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
 952     Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
 953     Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
 954     Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
 955
 956     evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
 957     evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
 958     evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
 959     evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
 960
 961     //loop track and get pid
 962     for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
 963       AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
 964
 965       if (!aodTrack){
 966         aodTrack->Delete();
 967         continue;
 968       }
 969
 970       Bool_t trkFlag = aodTrack->TestFilterBit(1); // TPC only tracks
 971       if(fFillDCA) trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
 972
 973       if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < 70) || !trkFlag){
 974         continue;
 975       }
 976
 977       Double_t b[2] = {-99., -99.};
 978       Double_t bCov[3] = {-99., -99., -99.};
 979       if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
 980         continue;
 981
 982       if (!fFillDCA && ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4)))
 983         continue;
 984
 985       if(fFillDCA && TMath::Abs(b[0]) > 3.0 && TMath::Abs(b[1]) > 3)
 986         continue;
 987
 988       // re-map the container in an array to do the analysis for V0A and V0C within a loop
 989       Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
 990       AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
 991       AliFlowVZEROQA *QA[2] = {fQA,fQA2};
 992
 993       Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
 994       AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
 995       AliFlowVZEROQA *QAv3[2] = {fQAv3,fQA2v3};
 996
 997       // Fill MC results
 998       if(fIsMC && mcArray){
 999         fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1000         Float_t tofMismProbMC = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1001
1002         Float_t xMC[5] = {iC,aodTrack->Charge(),1,evplaneMC,fPID->GetCurrentMask(1)&&tofMismProbMC < 0.5}; // to fill analysis v2 container
1003
1004         Float_t v2mc = TMath::Cos(2*(aodTrack->Phi() - evplaneMC));
1005
1006         fContAllChargesMC->Fill(0,aodTrack->Pt(),v2mc,xMC);
1007
1008         Int_t iS = TMath::Abs(((AliAODMCParticle*)mcArray->At(TMath::Abs(aodTrack->GetLabel())))->GetPdgCode());
1009         if(iS==11){
1010           fContAllChargesMC->Fill(4,aodTrack->Pt(),v2mc,xMC);
1011         }
1012         else if(iS==13){
1013           fContAllChargesMC->Fill(5,aodTrack->Pt(),v2mc,xMC);
1014         }
1015         else if(iS==211){
1016           fContAllChargesMC->Fill(1,aodTrack->Pt(),v2mc,xMC);
1017         }
1018         else if(iS==321){
1019           fContAllChargesMC->Fill(2,aodTrack->Pt(),v2mc,xMC);
1020         }
1021         else if(iS==2212){
1022           fContAllChargesMC->Fill(3,aodTrack->Pt(),v2mc,xMC);
1023         }
1024       }
1025
1026       for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1027
1028         if(fModulationDEDx) fPID->SetPsiCorrectionDeDx(evPlAngV0[iV0],evPlRes[iV0*8+iC]); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed)
1029
1030         Float_t v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1031         Float_t v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1032
1033         fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1034         Float_t dedx = fPID->GetDeDx();//aodTrack->GetTPCsignal();
1035         Float_t *probRead = fPID->GetProb();
1036         Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1037         Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1038         Float_t x[6] = {iC,aodTrack->Charge(),1,evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v2 container
1039         Float_t x3[6] = {iC,aodTrack->Charge(),1,evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v3 container
1040
1041         // in case fill DCA info
1042         if(fFillDCA){
1043           if(TMath::Abs(b[0]) > 0.1){
1044             x[5] = 1;
1045             x3[5] = 1;
1046           }
1047           if(TMath::Abs(b[0]) > 0.3){
1048             x[5] = 2;
1049             x3[5] = 2;
1050           }
1051           if(fIsMC && mcArray){
1052             if(!((AliAODMCParticle*)mcArray->At(TMath::Abs(aodTrack->GetLabel())))->IsPhysicalPrimary()){
1053               x[5] += 3;
1054               x3[5] += 3;
1055             }
1056           }
1057         }
1058
1059         // Fill no PID
1060         if(fV2) contV0[iV0]->Fill(0,aodTrack->Pt(),v2V0,x);
1061         if(fV3) contV0v3[iV0]->Fill(0,aodTrack->Pt(),v3V0,x3);
1062
1063
1064         Double_t dedxExp[8];
1065         Float_t tof = -1;
1066         Double_t inttimes[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
1067         Double_t expTOFsigma[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
1068
1069         Float_t nsigmaTPC[8];
1070         Float_t nsigmaTOF[8];
1071
1072         if(aodTrack->GetDetPid()){ // check the PID object is available
1073           for(Int_t iS=0;iS < 8;iS++){
1074             dedxExp[iS] = fPID->GetExpDeDx(aodTrack,iS);
1075             nsigmaTPC[iS] = (dedx - dedxExp[iS])/(dedxExp[iS]*0.07);
1076             //      printf("TPC %i = %f (%f %f)\n",iS, nsigmaTPC[iS],dedx, dedxExp[iS]);
1077           }
1078
1079           if(fPID->GetCurrentMask(1)){ // if TOF is present
1080             Float_t ptrack = aodTrack->P();
1081             tof = aodTrack->GetTOFsignal() - fPID->GetESDpid()->GetTOFResponse().GetStartTime(ptrack);
1082             aodTrack->GetIntegratedTimes(inttimes);
1083
1084             for(Int_t iS=5;iS < 8;iS++) // extra info for light nuclei
1085               inttimes[iS] = inttimes[0] / ptrack * mass[iS] * TMath::Sqrt(1+ptrack*ptrack/mass[iS]/mass[iS]);
1086
1087             for(Int_t iS=0;iS<8;iS++){
1088               expTOFsigma[iS] = fPID->GetESDpid()->GetTOFResponse().GetExpectedSigma(ptrack, inttimes[iS], mass[iS]);
1089               nsigmaTOF[iS] = (tof - inttimes[iS])/expTOFsigma[iS];
1090               //              printf("TOF %i = %f\n",iS, nsigmaTOF[iS]);
1091             }
1092           }
1093         }
1094
1095         Float_t deltaPhiV0 = aodTrack->Phi() - evPlAngV0[iV0];
1096         if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi();
1097         else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi();
1098         if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi();
1099         else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi();
1100
1101         Float_t deltaPhiV0v3 = aodTrack->Phi() - evPlAngV0v3[iV0];
1102         if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi();
1103         else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi();
1104         if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi();
1105         else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi();
1106
1107         // variable to fill QA container
1108         Float_t xQA[5] = {iC,aodTrack->Pt(), 0.0,deltaPhiV0,x[4]}; // v2
1109         Float_t xQA3[5] = {iC,aodTrack->Pt(), 0.0,deltaPhiV0v3,x[4]}; // v3
1110
1111         // extra QA TProfiles
1112         if(iV0==1 && aodTrack->Pt() < 20 && fPID->GetCurrentMask(0) && fPID->GetCurrentMask(1)){
1113           Float_t xQApid[2] = {iC,aodTrack->Pt()};
1114           fContQApid->Fill(0,nsigmaTPC[2],v2V0,xQApid); // v2 TPC (V0C) w.r.t pions
1115           fContQApid->Fill(1,nsigmaTOF[2],v2V0,xQApid); // v2 TOF (V0C) w.r.t. pions
1116           fContQApid->Fill(2,nsigmaTPC[3],v2V0,xQApid); // v2 TPC (V0C) w.r.t kaons
1117           fContQApid->Fill(3,nsigmaTOF[3],v2V0,xQApid); // v2 TOF (V0C) w.r.t. kaons
1118           fContQApid->Fill(4,nsigmaTPC[4],v2V0,xQApid); // v2 TPC (V0C) w.r.t protons
1119           fContQApid->Fill(5,nsigmaTOF[4],v2V0,xQApid); // v2 TOF (V0C) w.r.t. protons
1120         }
1121
1122         // QA fill
1123         if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid() || dedx < 10. || aodTrack->Pt() < 0 || aodTrack->Pt() > 7){}
1124         else{
1125           if(TMath::Abs(nsigmaTPC[2])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[2])<5))){ //pi
1126             xQA[2] = prob[2];
1127             xQA3[2] = xQA[2];
1128             if(fV2) QA[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA);
1129             if(fV3) QAv3[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA3);
1130           }
1131           if(TMath::Abs(nsigmaTPC[3])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[3])<5))){ //K
1132             xQA[2] = prob[3];
1133             xQA3[2] = xQA[2];
1134             if(fV2) QA[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA);
1135 //          if(fV3) QAv3[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA3);
1136           }
1137           if(TMath::Abs(nsigmaTPC[4])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[4])<5))){//p
1138             xQA[2] = prob[4];
1139             xQA3[2] = xQA[2];
1140             if(fV2) QA[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA);
1141 //          if(fV3) QAv3[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA3);
1142           }
1143           if(TMath::Abs(nsigmaTPC[0])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[0])<5))){//e
1144             xQA[2] = prob[0];
1145             xQA3[2] = xQA[2];
1146 //          if(fV2) QA[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA);
1147 //          if(fV3) QAv3[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA3);
1148           }
1149           if(TMath::Abs(nsigmaTPC[5])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[5])<5))){//d
1150             xQA[2] = prob[5];
1151             xQA3[2] = xQA[2];
1152             //    if(fV2) QA[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA);
1153             //    if(fV3) QAv3[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA3);
1154           }
1155           if(TMath::Abs(nsigmaTPC[6])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[6])<5))){//t
1156             xQA[2] = prob[6];
1157             xQA3[2] = xQA[2];
1158             //    if(fV2) QA[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA);
1159             //    if(fV3) QAv3[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA3);
1160           }
1161           if(TMath::Abs(nsigmaTPC[7])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[7])<5))){//He3
1162             xQA[2] = prob[7];
1163             xQA3[2] = xQA[2];
1164             //    if(fV2) QA[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA);
1165             //    if(fV3) QAv3[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA3);
1166           }
1167         }
1168
1169         //pid selection
1170         if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID and PID object strictly required (very important!!!!)
1171         else if(prob[2] > 0.6){ // pi
1172           x[2] = prob[2];
1173           x3[2] = x[2];
1174           if(TMath::Abs(nsigmaTPC[2]) < 5){ // TPC 5 sigma extra cut to accept the track
1175             if(fV2) contV0[iV0]->Fill(1,aodTrack->Pt(),v2V0,x);
1176             if(fV3) contV0v3[iV0]->Fill(1,aodTrack->Pt(),v3V0,x3);
1177             if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][0]->Fill(aodTrack->Pt(),b[0]);
1178             else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][0]->Fill(aodTrack->Pt(),b[0]);
1179           }
1180         }
1181         else if(prob[3] > 0.6){ // K
1182           x[2] = prob[3];
1183           x3[2] = x[2];
1184           if(TMath::Abs(nsigmaTPC[3]) < 5){
1185             if(fV2) contV0[iV0]->Fill(2,aodTrack->Pt(),v2V0,x);
1186             if(fV3) contV0v3[iV0]->Fill(2,aodTrack->Pt(),v3V0,x3);
1187             if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][1]->Fill(aodTrack->Pt(),b[0]);
1188             else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][1]->Fill(aodTrack->Pt(),b[0]);
1189           }
1190         }
1191         else if(prob[4] > 0.6){ // p
1192           x[2] = prob[4];
1193           x3[2] = x[2];
1194           if(TMath::Abs(nsigmaTPC[4]) < 5){
1195             if(fV2) contV0[iV0]->Fill(3,aodTrack->Pt(),v2V0,x);
1196             if(fV3) contV0v3[iV0]->Fill(3,aodTrack->Pt(),v3V0,x3);
1197             if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][2]->Fill(aodTrack->Pt(),b[0]);
1198             else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][2]->Fill(aodTrack->Pt(),b[0]);
1199           }
1200         }
1201         else if(prob[0] > 0.6){ // e
1202           x[2] = prob[0];
1203           x3[2] = x[2];
1204           if(TMath::Abs(nsigmaTPC[0]) < 5){
1205             if(fV2) contV0[iV0]->Fill(4,aodTrack->Pt(),v2V0,x);
1206             if(fV3) contV0v3[iV0]->Fill(4,aodTrack->Pt(),v3V0,x3);
1207             if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][3]->Fill(aodTrack->Pt(),b[0]);
1208             else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][3]->Fill(aodTrack->Pt(),b[0]);
1209           }
1210         }
1211         else if(prob[1] > 0.6){ // mu
1212           x[2] = prob[1];
1213           x3[2] = x[2];
1214           if(TMath::Abs(nsigmaTPC[1]) < 5){
1215             if(fV2) contV0[iV0]->Fill(8,aodTrack->Pt(),v2V0,x);
1216             if(fV3) contV0v3[iV0]->Fill(8,aodTrack->Pt(),v3V0,x3);
1217           }
1218         }
1219         else if(prob[5] > 0.6){ // d
1220           x[2] = prob[5];
1221           x3[2] = x[2];
1222           if(TMath::Abs(nsigmaTPC[5]) < 5){
1223             if(fV2) contV0[iV0]->Fill(5,aodTrack->Pt(),v2V0,x);
1224             if(fV3) contV0v3[iV0]->Fill(5,aodTrack->Pt(),v3V0,x3);
1225             if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][4]->Fill(aodTrack->Pt(),b[0]);
1226             else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][4]->Fill(aodTrack->Pt(),b[0]);
1227           }
1228         }
1229         else if(prob[6] > 0.6){ // t
1230           x[2] = prob[6];
1231           x3[2] = x[2];
1232           if(TMath::Abs(nsigmaTPC[6]) < 5){
1233             if(fV2) contV0[iV0]->Fill(6,aodTrack->Pt(),v2V0,x);
1234             if(fV3) contV0v3[iV0]->Fill(6,aodTrack->Pt(),v3V0,x3);
1235             if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][5]->Fill(aodTrack->Pt(),b[0]);
1236             else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][5]->Fill(aodTrack->Pt(),b[0]);
1237           }
1238         }
1239         else if(prob[7] > 0.6){ // He3
1240           x[2] = prob[7];
1241           x3[2] = x[2];
1242           if(TMath::Abs(nsigmaTPC[7]) < 5){
1243             if(fV2) contV0[iV0]->Fill(7,aodTrack->Pt()*2,v2V0,x);
1244             if(fV3) contV0v3[iV0]->Fill(7,aodTrack->Pt()*2,v3V0,x3);
1245             if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][6]->Fill(aodTrack->Pt(),b[0]);
1246             else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][6]->Fill(aodTrack->Pt(),b[0]);
1247           }
1248         }
1249
1250         if(x[4] > 0.5){ // if TOF was present redo TPC stand alone PID to check the PID in the same acceptance (PID mask = 2)
1251           fPID->ResetDetOR(1); // exclude TOF from PID
1252           tofMismProb = 0;
1253
1254           fPID->ComputeProb(aodTrack,fOutputAOD);
1255           dedx = fPID->GetDeDx();//aodTrack->GetTPCsignal();
1256           probRead = fPID->GetProb();
1257
1258           fPID->SetDetOR(1); // include TOF for PID
1259         }
1260         Float_t probTPC[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]}; // TPC stand alone prbabilities
1261
1262         //pid selection TPC S.A. with TOF matching
1263         x[4]*=2; // set the mask to 2 id TOF is present
1264         if(x[4]<1 || !(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID S.A. PID in TOF acceptance
1265         else if(probTPC[2] > 0.6){ // pi
1266           x[2] = probTPC[2];
1267           x3[2] = x[2];
1268           if(TMath::Abs(nsigmaTPC[2]) < 5){
1269             if(fV2) contV0[iV0]->Fill(1,aodTrack->Pt(),v2V0,x);
1270             if(fV3) contV0v3[iV0]->Fill(1,aodTrack->Pt(),v3V0,x3);
1271           }
1272         }
1273         else if(probTPC[3] > 0.6){ // K
1274           x[2] = probTPC[3];
1275           x3[2] = x[2];
1276           if(TMath::Abs(nsigmaTPC[3]) < 5){
1277             if(fV2) contV0[iV0]->Fill(2,aodTrack->Pt(),v2V0,x);
1278             if(fV3) contV0v3[iV0]->Fill(2,aodTrack->Pt(),v3V0,x3);
1279           }
1280         }
1281         else if(probTPC[4] > 0.6){ // p
1282           x[2] = probTPC[4];
1283           x3[2] = x[2];
1284           if(TMath::Abs(nsigmaTPC[4]) < 5){
1285             if(fV2) contV0[iV0]->Fill(3,aodTrack->Pt(),v2V0,x);
1286             if(fV3) contV0v3[iV0]->Fill(3,aodTrack->Pt(),v3V0,x3);
1287           }
1288         }
1289         else if(probTPC[0] > 0.6){ // e
1290           x[2] = probTPC[0];
1291           x3[2] = x[2];
1292           if(TMath::Abs(nsigmaTPC[0]) < 5){
1293             if(fV2) contV0[iV0]->Fill(4,aodTrack->Pt(),v2V0,x);
1294             if(fV3) contV0v3[iV0]->Fill(4,aodTrack->Pt(),v3V0,x3);
1295           }
1296         }
1297         else if(probTPC[1] > 0.6){ // mu
1298           x[2] = probTPC[1];
1299           x3[2] = x[2];
1300           if(TMath::Abs(nsigmaTPC[1]) < 5){
1301             if(fV2) contV0[iV0]->Fill(8,aodTrack->Pt(),v2V0,x);
1302             if(fV3) contV0v3[iV0]->Fill(8,aodTrack->Pt(),v3V0,x3);
1303           }
1304         }
1305         else if(probTPC[5] > 0.6){ // d
1306           x[2] = probTPC[5];
1307           x3[2] = x[2];
1308           if(TMath::Abs(nsigmaTPC[5]) < 5){
1309             if(fV2) contV0[iV0]->Fill(5,aodTrack->Pt(),v2V0,x);
1310             if(fV3) contV0v3[iV0]->Fill(5,aodTrack->Pt(),v3V0,x3);
1311           }
1312         }
1313         else if(probTPC[6] > 0.6){ // t
1314           x[2] = probTPC[6];
1315           x3[2] = x[2];
1316           if(TMath::Abs(nsigmaTPC[6]) < 5){
1317             if(fV2) contV0[iV0]->Fill(6,aodTrack->Pt(),v2V0,x);
1318             if(fV3) contV0v3[iV0]->Fill(6,aodTrack->Pt(),v3V0,x3);
1319           }
1320         }
1321         else if(probTPC[7] > 0.6){ // He3
1322           x[2] = probTPC[7];
1323           x3[2] = x[2];
1324           if(TMath::Abs(nsigmaTPC[7]) < 5){
1325             if(fV2) contV0[iV0]->Fill(7,aodTrack->Pt()*2,v2V0,x);
1326             if(fV3) contV0v3[iV0]->Fill(7,aodTrack->Pt()*2,v3V0,x3);
1327           }
1328         }
1329       } // end side loop
1330     } // end track loop
1331
1332     // Fill EP distribution histograms
1333     if(fV2) fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
1334     if(fV2) fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
1335
1336     if(fV3) fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
1337     if(fV3) fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
1338
1339     // TPC EP needed for resolution studies (TPC subevent)
1340     Double_t Qx2 = 0, Qy2 = 0;
1341     Double_t Qx3 = 0, Qy3 = 0;
1342
1343     for(Int_t iT = 0; iT < nAODTracks; iT++) {
1344
1345       AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
1346
1347       if (!aodTrack){
1348         aodTrack->Delete();
1349         continue;
1350       }
1351
1352       Bool_t trkFlag = aodTrack->TestFilterBit(1);
1353
1354       if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70)  || !trkFlag)
1355         continue;
1356
1357       Double_t b[2] = {-99., -99.};
1358       Double_t bCov[3] = {-99., -99., -99.};
1359       if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
1360         continue;
1361
1362       if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
1363         continue;
1364
1365       Qx2 += TMath::Cos(2*aodTrack->Phi());
1366       Qy2 += TMath::Sin(2*aodTrack->Phi());
1367       Qx3 += TMath::Cos(3*aodTrack->Phi());
1368       Qy3 += TMath::Sin(3*aodTrack->Phi());
1369
1370     }
1371
1372     evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
1373     evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
1374
1375     // Fill histograms needed for resolution evaluation
1376     if(fV2) fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
1377     if(fV2) fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
1378     if(fV2) fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
1379
1380     if(fV3) fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
1381     if(fV3) fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
1382     if(fV3) fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
1383   }
1384
1385 }
1386
1387 //_____________________________________________________________________________
1388 Float_t AliAnalysisTaskVnV0::GetVertex(AliAODEvent* aod) const
1389 {
1390
1391   Float_t zvtx = -999;
1392
1393   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
1394   if (!vtxAOD)
1395     return zvtx;
1396   if(vtxAOD->GetNContributors()>0)
1397     zvtx = vtxAOD->GetZ();
1398
1399   return zvtx;
1400 }
1401 //_____________________________________________________________________________
1402 void AliAnalysisTaskVnV0::Terminate(Option_t *)
1403 {
1404   // Terminate loop
1405   Printf("Terminate()");
1406 }
1407 //_____________________________________________________________________________
1408 void AliAnalysisTaskVnV0::OpenInfoCalbration(Int_t run){
1409     TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
1410     TFile *foadb = TFile::Open(oadbfilename.Data());
1411
1412     if(!foadb){
1413         printf("OADB file %s cannot be opened\n",oadbfilename.Data());
1414         return;
1415     }
1416
1417     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1418     if(!cont){
1419         printf("OADB object hMultV0BefCorr is not available in the file\n");
1420         return;
1421     }
1422
1423     if(!(cont->GetObject(run))){
1424         printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
1425         run = 137366;
1426     }
1427     fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
1428
1429     TF1 *fpol0 = new TF1("fpol0","pol0");
1430     fMultV0->Fit(fpol0,"","",0,31);
1431     fV0Cpol = fpol0->GetParameter(0);
1432     fMultV0->Fit(fpol0,"","",32,64);
1433     fV0Apol = fpol0->GetParameter(0);
1434
1435     for(Int_t iside=0;iside<2;iside++){
1436         for(Int_t icoord=0;icoord<2;icoord++){
1437             for(Int_t i=0;i  < nCentrBin;i++){
1438                 char namecont[100];
1439                 if(iside==0 && icoord==0)
1440                   snprintf(namecont,100,"hQxc2_%i",i);
1441                 else if(iside==1 && icoord==0)
1442                   snprintf(namecont,100,"hQxa2_%i",i);
1443                 else if(iside==0 && icoord==1)
1444                   snprintf(namecont,100,"hQyc2_%i",i);
1445                 else if(iside==1 && icoord==1)
1446                   snprintf(namecont,100,"hQya2_%i",i);
1447
1448                 cont = (AliOADBContainer*) foadb->Get(namecont);
1449                 if(!cont){
1450                     printf("OADB object %s is not available in the file\n",namecont);
1451                     return;
1452                 }
1453
1454                 if(!(cont->GetObject(run))){
1455                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1456                     run = 137366;
1457                 }
1458                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1459                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1460
1461                 //for v3
1462                 if(iside==0 && icoord==0)
1463                   snprintf(namecont,100,"hQxc3_%i",i);
1464                 else if(iside==1 && icoord==0)
1465                   snprintf(namecont,100,"hQxa3_%i",i);
1466                 else if(iside==0 && icoord==1)
1467                   snprintf(namecont,100,"hQyc3_%i",i);
1468                 else if(iside==1 && icoord==1)
1469                   snprintf(namecont,100,"hQya3_%i",i);
1470
1471                 cont = (AliOADBContainer*) foadb->Get(namecont);
1472                 if(!cont){
1473                     printf("OADB object %s is not available in the file\n",namecont);
1474                     return;
1475                 }
1476
1477                 if(!(cont->GetObject(run))){
1478                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1479                     run = 137366;
1480                 }
1481                 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1482                 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1483
1484             }
1485         }
1486     }
1487 }