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 #include "AliESDtrackCuts.h"
  24 #include "AliESDVertex.h"
  25 #include "AliEventplane.h"
  26 #include "TProfile2D.h"
  27
  28 // STL includes
  29 //#include <iostream>
  30 //using namespace std;
  31
  32 ClassImp(AliAnalysisTaskVnV0)
  33 Bool_t AliAnalysisTaskVnV0::fgIsPsiComputed = kFALSE;
  34 Float_t AliAnalysisTaskVnV0::fgPsi2v0a=999.;
  35 Float_t AliAnalysisTaskVnV0::fgPsi2v0c=999.;
  36 Float_t AliAnalysisTaskVnV0::fgPsi2tpc=999.;
  37 Float_t AliAnalysisTaskVnV0::fgPsi3v0a=999.;
  38 Float_t AliAnalysisTaskVnV0::fgPsi3v0c=999.;
  39 Float_t AliAnalysisTaskVnV0::fgPsi3tpc=999.;
  40 Float_t AliAnalysisTaskVnV0::fgPsi2v0aMC=999.;
  41 Float_t AliAnalysisTaskVnV0::fgPsi2v0cMC=999.;
  42 Float_t AliAnalysisTaskVnV0::fgPsi2tpcMC=999.;
  43 Float_t AliAnalysisTaskVnV0::fgPsi3v0aMC=999.;
  44 Float_t AliAnalysisTaskVnV0::fgPsi3v0cMC=999.;
  45 Float_t AliAnalysisTaskVnV0::fgPsi3tpcMC=999.;
  46
  47 //_____________________________________________________________________________
  48 AliAnalysisTaskVnV0::AliAnalysisTaskVnV0():
  49   AliAnalysisTaskSE(),
  50   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
  51   fEtaCut(0.8),   // cut on |eta| < fEtaCut
  52   fMinPt(0.15),   // cut on pt > fMinPt
  53   fMinDistV0(0),
  54   fMaxDistV0(100),
  55   fV2(kTRUE),
  56   fV3(kTRUE),
  57   fIsMC(kFALSE),
  58   fQAsw(kFALSE),
  59   fIsAfter2011(kFALSE),
  60   fRun(-1),
  61   fNcluster(70),
  62   fList(new TList()),
  63   fList2(new TList()),
  64   fList3(new TList()),
  65   fList4(new TList()),
  66   fMultV0(NULL),
  67   fV0Cpol(100),
  68   fV0Apol(100),
  69   fHResTPCv0A2(NULL),
  70   fHResTPCv0C2(NULL),
  71   fHResv0Cv0A2(NULL),
  72   fHResTPCv0A3(NULL),
  73   fHResTPCv0C3(NULL),
  74   fHResv0Cv0A3(NULL),
  75   fPhiRPv0A(NULL),
  76   fPhiRPv0C(NULL),
  77   fPhiRPv0Av3(NULL),
  78   fPhiRPv0Cv3(NULL),
  79   fQA(NULL),
  80   fQA2(NULL),
  81   fQAv3(NULL),
  82   fQA2v3(NULL),
  83   fPID(new AliFlowBayesianPID()),
  84   fTree(NULL),
  85   fCentrality(-1),
  86   evPlAngV0ACor2(0),
  87   evPlAngV0CCor2(0),
  88   evPlAng2(0),
  89   evPlAngV0ACor3(0),
  90   evPlAngV0CCor3(0),
  91   evPlAng3(0),
  92   fContAllChargesV0A(NULL),
  93   fContAllChargesV0C(NULL),
  94   fContAllChargesV0Av3(NULL),
  95   fContAllChargesV0Cv3(NULL),
  96   fContAllChargesMC(NULL),
  97   fHResMA2(NULL),
  98   fHResMC2(NULL),
  99   fHResAC2(NULL),
 100   fHResMA3(NULL),
 101   fHResMC3(NULL),
 102   fHResAC3(NULL),
 103   fContAllChargesMCA(NULL),
 104   fContAllChargesMCC(NULL),
 105   fContAllChargesMCAv3(NULL),
 106   fContAllChargesMCCv3(NULL),
 107   fFillDCA(kFALSE),
 108   fContQApid(NULL),
 109   fModulationDEDx(kFALSE),
 110   fZvtx(0.),
 111   fNK0s(0),
 112   fNpiPos(0),
 113   fNpiNeg(0),
 114   fHKsPhi(NULL),
 115   fHKsPhiEP(NULL),
 116   fHK0sMass(NULL),
 117   fHK0sMass2(NULL),
 118   fHK0vsLambda(NULL),
 119   fHctauPtEP(NULL),
 120   fHctauAt1EP(NULL),
 121   fCutsDaughter(NULL)
 122 {
 123   // Default constructor (should not be used)
 124   fList->SetName("resultsV2");
 125   fList2->SetName("resultsV3");
 126   fList3->SetName("resultsMC");
 127   fList4->SetName("QA");
 128
 129   fList->SetOwner(kTRUE);
 130   fList2->SetOwner(kTRUE);
 131   fList3->SetOwner(kTRUE);
 132   fList4->SetOwner(kTRUE);
 133
 134   fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h
 135 }
 136
 137 //______________________________________________________________________________
 138 AliAnalysisTaskVnV0::AliAnalysisTaskVnV0(const char *name):
 139   AliAnalysisTaskSE(name),
 140   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
 141   fEtaCut(0.8),   // cut on |eta| < fEtaCut
 142   fMinPt(0.15),   // cut on pt > fMinPt
 143   fMinDistV0(0),
 144   fMaxDistV0(100),
 145   fV2(kTRUE),
 146   fV3(kTRUE),
 147   fIsMC(kFALSE),
 148   fQAsw(kFALSE),
 149   fIsAfter2011(kFALSE),
 150   fRun(-1),
 151   fNcluster(70),
 152   fList(new TList()),
 153   fList2(new TList()),
 154   fList3(new TList()),
 155   fList4(new TList()),
 156   fMultV0(NULL),
 157   fV0Cpol(100),
 158   fV0Apol(100),
 159   fHResTPCv0A2(NULL),
 160   fHResTPCv0C2(NULL),
 161   fHResv0Cv0A2(NULL),
 162   fHResTPCv0A3(NULL),
 163   fHResTPCv0C3(NULL),
 164   fHResv0Cv0A3(NULL),
 165   fPhiRPv0A(NULL),
 166   fPhiRPv0C(NULL),
 167   fPhiRPv0Av3(NULL),
 168   fPhiRPv0Cv3(NULL),
 169   fQA(NULL),
 170   fQA2(NULL),
 171   fQAv3(NULL),
 172   fQA2v3(NULL),
 173   fPID(new AliFlowBayesianPID()),
 174   fTree(NULL),
 175   fCentrality(-1),
 176   evPlAngV0ACor2(0),
 177   evPlAngV0CCor2(0),
 178   evPlAng2(0),
 179   evPlAngV0ACor3(0),
 180   evPlAngV0CCor3(0),
 181   evPlAng3(0),
 182   fContAllChargesV0A(NULL),
 183   fContAllChargesV0C(NULL),
 184   fContAllChargesV0Av3(NULL),
 185   fContAllChargesV0Cv3(NULL),
 186   fContAllChargesMC(NULL),
 187   fHResMA2(NULL),
 188   fHResMC2(NULL),
 189   fHResAC2(NULL),
 190   fHResMA3(NULL),
 191   fHResMC3(NULL),
 192   fHResAC3(NULL),
 193   fContAllChargesMCA(NULL),
 194   fContAllChargesMCC(NULL),
 195   fContAllChargesMCAv3(NULL),
 196   fContAllChargesMCCv3(NULL),
 197   fFillDCA(kFALSE),
 198   fContQApid(NULL),
 199   fModulationDEDx(kFALSE),
 200   fZvtx(0.),
 201   fNK0s(0),
 202   fNpiPos(0),
 203   fNpiNeg(0),
 204   fHKsPhi(NULL),
 205   fHKsPhiEP(NULL),
 206   fHK0sMass(NULL),
 207   fHK0sMass2(NULL),
 208   fHK0vsLambda(NULL),
 209   fHctauPtEP(NULL),
 210   fHctauAt1EP(NULL),
 211   fCutsDaughter(NULL)
 212 {
 213
 214   DefineOutput(1, TList::Class());
 215   DefineOutput(2, TList::Class());
 216   DefineOutput(3, TList::Class());
 217   DefineOutput(4, TList::Class());
 218
 219   // Output slot #1 writes into a TTree
 220   fList->SetName("resultsV2");
 221   fList2->SetName("resultsV3");
 222   fList3->SetName("resultsMC");
 223   fList4->SetName("QA");
 224
 225   fList->SetOwner(kTRUE);
 226   fList2->SetOwner(kTRUE);
 227   fList3->SetOwner(kTRUE);
 228   fList4->SetOwner(kTRUE);
 229
 230   fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h
 231 }
 232 //_____________________________________________________________________________
 233 AliAnalysisTaskVnV0::~AliAnalysisTaskVnV0()
 234 {
 235
 236 }
 237
 238 //______________________________________________________________________________
 239 void AliAnalysisTaskVnV0::UserCreateOutputObjects()
 240 {
 241
 242   if(fIsMC) fPID->SetMC(kTRUE);
 243
 244
 245   // Tree for EP debug (comment the adding to v2 list id not needed)
 246   fTree = new TTree("tree","tree");
 247   fTree->Branch("evPlAngV0ACor2",&evPlAngV0ACor2,"evPlAngV0ACor2/F");
 248   fTree->Branch("evPlAngV0CCor2",&evPlAngV0CCor2,"evPlAngV0CCor2/F");
 249   fTree->Branch("evPlAng2",&evPlAng2,"evPlAng2/F");
 250   fTree->Branch("fCentrality",&fCentrality,"fCentrality/F");
 251   fTree->Branch("evPlAngV0ACor3",&evPlAngV0ACor3,"evPlAngV0ACor3/F");
 252   fTree->Branch("evPlAngV0CCor3",&evPlAngV0CCor3,"evPlAngV0CCor3/F");
 253   fTree->Branch("evPlAng3",&evPlAng3,"evPlAng3/F");
 254
 255
 256   // Container analyses (different steps mean different species)
 257   const Int_t nPtBinsTOF = 45;
 258   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};
 259   const Int_t nCentrTOF = nCentrBin;
 260   const Int_t nPsiTOF = 10;
 261   const Int_t nChargeBinsTOFres = 2;
 262   const Int_t nCentrTOFres = nCentrBin;
 263   const Int_t nProbTOFres = 4;
 264   const Int_t nPsiTOFres = 10;
 265   const Int_t nMaskPID = 3;
 266
 267   Int_t nDCABin = 1; // put to 1 not to store this info
 268   if(fFillDCA) nDCABin = 3;
 269   if(fIsMC && nDCABin>1)  nDCABin = 6;
 270   /*
 271     0 = DCAxy < 2.4 && all (or Physical primary if MC)
 272     1 = DCAxy > 2.4 && all (or Physical primary if MC)
 273     2 = DCAxy < 2.4 && not Physical Primary for MC
 274     3 = DCAxy > 2.4 && not Physical Primary for MC
 275   */
 276
 277   Int_t binsTOF[6] = {nCentrTOFres,nChargeBinsTOFres,nProbTOFres,nPsiTOFres,nMaskPID,nDCABin};
 278   Int_t binsTOFmc[5] = {nCentrTOFres,nChargeBinsTOFres,1,nPsiTOFres,2};
 279   Int_t binsTOFmcPureMC[5] = {nCentrTOFres,nChargeBinsTOFres,1,nPsiTOFres,1};
 280
 281   // v2 container
 282   fContAllChargesV0A = new AliFlowVZEROResults("v2A",6,binsTOF);
 283   fContAllChargesV0A->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 284   fContAllChargesV0A->SetVarRange(1,-1.5,1.5);  // charge
 285   fContAllChargesV0A->SetVarRange(2,0.6,1.0001);// prob
 286   fContAllChargesV0A->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
 287   fContAllChargesV0A->SetVarRange(4,-0.5,2.5); // pid mask
 288   fContAllChargesV0A->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
 289   fContAllChargesV0A->SetVarName(0,"centrality");
 290   fContAllChargesV0A->SetVarName(1,"charge");
 291   fContAllChargesV0A->SetVarName(2,"prob");
 292   fContAllChargesV0A->SetVarName(3,"#Psi");
 293   fContAllChargesV0A->SetVarName(4,"PIDmask");
 294   fContAllChargesV0A->SetVarName(5,"DCAbin");
 295   if(fV2) fContAllChargesV0A->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 296   if(fV2) fContAllChargesV0A->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 297   if(fV2) fContAllChargesV0A->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 298   if(fV2) fContAllChargesV0A->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 299   if(fV2) fContAllChargesV0A->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 300   if(fV2) fContAllChargesV0A->AddSpecies("d",nPtBinsTOF,binsPtTOF);
 301   if(fV2) fContAllChargesV0A->AddSpecies("t",nPtBinsTOF,binsPtTOF);
 302   if(fV2) fContAllChargesV0A->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
 303   if(fV2) fContAllChargesV0A->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 304   if(fV2) fContAllChargesV0A->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
 305   if(fV2) fContAllChargesV0A->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
 306   if(fV2) fContAllChargesV0A->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
 307   if(fV2) fContAllChargesV0A->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
 308
 309   fContAllChargesV0C = new AliFlowVZEROResults("v2C",6,binsTOF);
 310   fContAllChargesV0C->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 311   fContAllChargesV0C->SetVarRange(1,-1.5,1.5);  // charge
 312   fContAllChargesV0C->SetVarRange(2,0.6,1.0001);// prob
 313   fContAllChargesV0C->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
 314   fContAllChargesV0C->SetVarRange(4,-0.5,2.5); // pid mask
 315   fContAllChargesV0C->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
 316   fContAllChargesV0C->SetVarName(0,"centrality");
 317   fContAllChargesV0C->SetVarName(1,"charge");
 318   fContAllChargesV0C->SetVarName(2,"prob");
 319   fContAllChargesV0C->SetVarName(3,"#Psi");
 320   fContAllChargesV0C->SetVarName(4,"PIDmask");
 321   fContAllChargesV0C->SetVarName(5,"DCAbin");
 322   if(fV2) fContAllChargesV0C->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 323   if(fV2) fContAllChargesV0C->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 324   if(fV2) fContAllChargesV0C->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 325   if(fV2) fContAllChargesV0C->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 326   if(fV2) fContAllChargesV0C->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 327   if(fV2) fContAllChargesV0C->AddSpecies("d",nPtBinsTOF,binsPtTOF);
 328   if(fV2) fContAllChargesV0C->AddSpecies("t",nPtBinsTOF,binsPtTOF);
 329   if(fV2) fContAllChargesV0C->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
 330   if(fV2) fContAllChargesV0C->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 331   if(fV2) fContAllChargesV0C->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
 332   if(fV2) fContAllChargesV0C->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
 333   if(fV2) fContAllChargesV0C->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
 334   if(fV2) fContAllChargesV0C->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
 335
 336   fList->Add(fContAllChargesV0A);
 337   fList->Add(fContAllChargesV0C);
 338
 339   fHctauPtEP = new TProfile2D("hctauPtEP","K^{0}_{s} decay length;p_{T} (GeV/#it{c});#Delta#phi (rad)",40,0,5,10,-TMath::Pi(),TMath::Pi());
 340   fHctauAt1EP = new TH2F("hctauAt1EP","K^{0}_{s} decay length at 1 GeV/#it{c};c#tau (cm);#Delta#phi (rad)",50,0,50,10,-TMath::Pi(),TMath::Pi());
 341   // added at the end
 342
 343   if(fIsMC && fV2){
 344     fContAllChargesMC = new AliFlowVZEROResults("v2mc",5,binsTOFmc);
 345     fContAllChargesMC->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 346     fContAllChargesMC->SetVarRange(1,-1.5,1.5);  // charge
 347     fContAllChargesMC->SetVarRange(2,0.6,1.0001);// prob
 348     fContAllChargesMC->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
 349     fContAllChargesMC->SetVarRange(4,-0.5,1.5); // pid mask
 350     fContAllChargesMC->SetVarName(0,"centrality");
 351     fContAllChargesMC->SetVarName(1,"charge");
 352     fContAllChargesMC->SetVarName(2,"prob");
 353     fContAllChargesMC->SetVarName(3,"#Psi");
 354     fContAllChargesMC->SetVarName(4,"PIDmask");
 355     fContAllChargesMC->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 356     fContAllChargesMC->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 357     fContAllChargesMC->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 358     fContAllChargesMC->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 359     fContAllChargesMC->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 360     fContAllChargesMC->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 361     fList3->Add(fContAllChargesMC);
 362
 363     fContAllChargesMCA = new AliFlowVZEROResults("v2mcA",5,binsTOFmcPureMC);
 364     fContAllChargesMCA->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 365     fContAllChargesMCA->SetVarRange(1,-1.5,1.5);  // charge
 366     fContAllChargesMCA->SetVarRange(2,0.6,1.0001);// prob
 367     fContAllChargesMCA->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
 368     fContAllChargesMCA->SetVarRange(4,-0.5,1.5); // pid mask
 369     fContAllChargesMCA->SetVarName(0,"centrality");
 370     fContAllChargesMCA->SetVarName(1,"charge");
 371     fContAllChargesMCA->SetVarName(2,"prob");
 372     fContAllChargesMCA->SetVarName(3,"#Psi");
 373     fContAllChargesMCA->SetVarName(4,"PIDmask");
 374     fContAllChargesMCA->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 375     fContAllChargesMCA->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 376     fContAllChargesMCA->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 377     fContAllChargesMCA->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 378     fContAllChargesMCA->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 379     fContAllChargesMCA->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 380     fList3->Add(fContAllChargesMCA);
 381
 382     fContAllChargesMCC = new AliFlowVZEROResults("v2mcC",5,binsTOFmcPureMC);
 383     fContAllChargesMCC->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 384     fContAllChargesMCC->SetVarRange(1,-1.5,1.5);  // charge
 385     fContAllChargesMCC->SetVarRange(2,0.6,1.0001);// prob
 386     fContAllChargesMCC->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
 387     fContAllChargesMCC->SetVarRange(4,-0.5,1.5); // pid mask
 388     fContAllChargesMCC->SetVarName(0,"centrality");
 389     fContAllChargesMCC->SetVarName(1,"charge");
 390     fContAllChargesMCC->SetVarName(2,"prob");
 391     fContAllChargesMCC->SetVarName(3,"#Psi");
 392     fContAllChargesMCC->SetVarName(4,"PIDmask");
 393     fContAllChargesMCC->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 394     fContAllChargesMCC->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 395     fContAllChargesMCC->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 396     fContAllChargesMCC->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 397     fContAllChargesMCC->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 398     fContAllChargesMCC->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 399     fList3->Add(fContAllChargesMCC);
 400   }
 401
 402   // v3 container
 403   fContAllChargesV0Av3 = new AliFlowVZEROResults("v3A",6,binsTOF);
 404   fContAllChargesV0Av3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 405   fContAllChargesV0Av3->SetVarRange(1,-1.5,1.5);  // charge
 406   fContAllChargesV0Av3->SetVarRange(2,0.6,1.0001);// prob
 407   fContAllChargesV0Av3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
 408   fContAllChargesV0Av3->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
 409   fContAllChargesV0Av3->SetVarRange(4,-0.5,2.5); // pid mask
 410   fContAllChargesV0Av3->SetVarName(0,"centrality");
 411   fContAllChargesV0Av3->SetVarName(1,"charge");
 412   fContAllChargesV0Av3->SetVarName(2,"prob");
 413   fContAllChargesV0Av3->SetVarName(3,"#Psi");
 414   fContAllChargesV0Av3->SetVarName(4,"PIDmask");
 415   fContAllChargesV0Av3->SetVarName(5,"DCAbin");
 416   if(fV3) fContAllChargesV0Av3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 417   if(fV3) fContAllChargesV0Av3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 418   if(fV3) fContAllChargesV0Av3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 419   if(fV3) fContAllChargesV0Av3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 420   if(fV3) fContAllChargesV0Av3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 421   if(fV3) fContAllChargesV0Av3->AddSpecies("d",nPtBinsTOF,binsPtTOF);
 422   if(fV3) fContAllChargesV0Av3->AddSpecies("t",nPtBinsTOF,binsPtTOF);
 423   if(fV3) fContAllChargesV0Av3->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
 424   if(fV3) fContAllChargesV0Av3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 425   if(fV3) fContAllChargesV0Av3->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
 426   if(fV3) fContAllChargesV0Av3->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
 427   if(fV3) fContAllChargesV0Av3->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
 428   if(fV3) fContAllChargesV0Av3->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
 429
 430   fContAllChargesV0Cv3 = new AliFlowVZEROResults("v3C",6,binsTOF);
 431   fContAllChargesV0Cv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 432   fContAllChargesV0Cv3->SetVarRange(1,-1.5,1.5);  // charge
 433   fContAllChargesV0Cv3->SetVarRange(2,0.6,1.0001);// prob
 434   fContAllChargesV0Cv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
 435   fContAllChargesV0Cv3->SetVarRange(4,-0.5,2.5); // pid mask
 436   fContAllChargesV0Cv3->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
 437   fContAllChargesV0Cv3->SetVarName(0,"centrality");
 438   fContAllChargesV0Cv3->SetVarName(1,"charge");
 439   fContAllChargesV0Cv3->SetVarName(2,"prob");
 440   fContAllChargesV0Cv3->SetVarName(3,"#Psi");
 441   fContAllChargesV0Cv3->SetVarName(4,"PIDmask");
 442   fContAllChargesV0Cv3->SetVarName(5,"DCAbin");
 443   if(fV3) fContAllChargesV0Cv3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 444   if(fV3) fContAllChargesV0Cv3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 445   if(fV3) fContAllChargesV0Cv3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 446   if(fV3) fContAllChargesV0Cv3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 447   if(fV3) fContAllChargesV0Cv3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 448   if(fV3) fContAllChargesV0Cv3->AddSpecies("d",nPtBinsTOF,binsPtTOF);
 449   if(fV3) fContAllChargesV0Cv3->AddSpecies("t",nPtBinsTOF,binsPtTOF);
 450   if(fV3) fContAllChargesV0Cv3->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
 451   if(fV3) fContAllChargesV0Cv3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 452   if(fV3) fContAllChargesV0Cv3->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
 453   if(fV3) fContAllChargesV0Cv3->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
 454   if(fV3) fContAllChargesV0Cv3->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
 455   if(fV3) fContAllChargesV0Cv3->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
 456
 457   fList2->Add(fContAllChargesV0Av3);
 458   fList2->Add(fContAllChargesV0Cv3);
 459
 460   if(fIsMC && fV3){
 461     fContAllChargesMCAv3 = new AliFlowVZEROResults("v3mcA",5,binsTOFmcPureMC);
 462     fContAllChargesMCAv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 463     fContAllChargesMCAv3->SetVarRange(1,-1.5,1.5);  // charge
 464     fContAllChargesMCAv3->SetVarRange(2,0.6,1.0001);// prob
 465     fContAllChargesMCAv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
 466     fContAllChargesMCAv3->SetVarRange(4,-0.5,1.5); // pid mask
 467     fContAllChargesMCAv3->SetVarName(0,"centrality");
 468     fContAllChargesMCAv3->SetVarName(1,"charge");
 469     fContAllChargesMCAv3->SetVarName(2,"prob");
 470     fContAllChargesMCAv3->SetVarName(3,"#Psi");
 471     fContAllChargesMCAv3->SetVarName(4,"PIDmask");
 472     fContAllChargesMCAv3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 473     fContAllChargesMCAv3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 474     fContAllChargesMCAv3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 475     fContAllChargesMCAv3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 476     fContAllChargesMCAv3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 477     fContAllChargesMCAv3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 478     fList3->Add(fContAllChargesMCAv3);
 479
 480     fContAllChargesMCCv3 = new AliFlowVZEROResults("v3mcC",5,binsTOFmcPureMC);
 481     fContAllChargesMCCv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 482     fContAllChargesMCCv3->SetVarRange(1,-1.5,1.5);  // charge
 483     fContAllChargesMCCv3->SetVarRange(2,0.6,1.0001);// prob
 484     fContAllChargesMCCv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
 485     fContAllChargesMCCv3->SetVarRange(4,-0.5,1.5); // pid mask
 486     fContAllChargesMCCv3->SetVarName(0,"centrality");
 487     fContAllChargesMCCv3->SetVarName(1,"charge");
 488     fContAllChargesMCCv3->SetVarName(2,"prob");
 489     fContAllChargesMCCv3->SetVarName(3,"#Psi");
 490     fContAllChargesMCCv3->SetVarName(4,"PIDmask");
 491     fContAllChargesMCCv3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
 492     fContAllChargesMCCv3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
 493     fContAllChargesMCCv3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
 494     fContAllChargesMCCv3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
 495     fContAllChargesMCCv3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
 496     fContAllChargesMCCv3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
 497     fList3->Add(fContAllChargesMCCv3);
 498   }
 499
 500   // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
 501   // v2
 502   fHResTPCv0A2 = new TProfile("hResTPCv0A2","",nCentrBin,0,nCentrBin);
 503   fHResTPCv0C2 = new TProfile("hResTPCv0C2","",nCentrBin,0,nCentrBin);
 504   fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",nCentrBin,0,nCentrBin);
 505
 506   fList->Add(fHResTPCv0A2);
 507   fList->Add(fHResTPCv0C2);
 508   fList->Add(fHResv0Cv0A2);
 509
 510   // v3
 511   fHResTPCv0A3 = new TProfile("hResTPCv0A3","",nCentrBin,0,nCentrBin);
 512   fHResTPCv0C3 = new TProfile("hResTPCv0C3","",nCentrBin,0,nCentrBin);
 513   fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",nCentrBin,0,nCentrBin);
 514
 515   fList2->Add(fHResTPCv0A3);
 516   fList2->Add(fHResTPCv0C3);
 517   fList2->Add(fHResv0Cv0A3);
 518
 519   // MC as in the dataEP resolution (but using MC tracks)
 520   if(fIsMC && fV3){
 521     fHResMA2 = new TProfile("hResMA2","",nCentrBin,0,nCentrBin);
 522     fHResMC2 = new TProfile("hResMC2","",nCentrBin,0,nCentrBin);
 523     fHResAC2 = new TProfile("hResAC2","",nCentrBin,0,nCentrBin);
 524     fList3->Add(fHResMA2);
 525     fList3->Add(fHResMC2);
 526     fList3->Add(fHResAC2);
 527   }
 528   if(fIsMC && fV3){
 529     fHResMA3 = new TProfile("hResMA3","",nCentrBin,0,nCentrBin);
 530     fHResMC3 = new TProfile("hResMC3","",nCentrBin,0,nCentrBin);
 531     fHResAC3 = new TProfile("hResAC3","",nCentrBin,0,nCentrBin);
 532     fList3->Add(fHResMA3);
 533     fList3->Add(fHResMC3);
 534     fList3->Add(fHResAC3);
 535   }
 536
 537
 538   // V0A and V0C event plane distributions
 539   //v2
 540   fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
 541   fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
 542
 543   //v3
 544   fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
 545   fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
 546
 547   // QA container
 548   // v2
 549   const Int_t nDETsignal = 50;
 550   Double_t binDETsignal[nDETsignal+1];
 551   for(Int_t i=0;i<nDETsignal+1;i++){
 552     binDETsignal[i] = -5 + i*10. / nDETsignal;
 553   }
 554 //   const Int_t nEta = 5;
 555 //   Double_t binEta[nEta+1];
 556 //   for(Int_t i=0;i<nEta+1;i++){
 557 //     binEta[i] = -1 + i*2. / nEta;
 558 //   }
 559
 560   const Int_t nDeltaPhi = 5;
 561   const Int_t nDeltaPhiV3 = 7;
 562
 563   Int_t binsQA[5] = {nCentrTOF,7,5,nDeltaPhi,2};
 564   Int_t binsQAv3[5] = {nCentrTOF,7,5,nDeltaPhiV3,2};
 565
 566
 567   fQA = new AliFlowVZEROQA("v2AQA",5,binsQA);
 568   fQA->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 569   fQA->SetVarRange(1,0,7);  // pt
 570   fQA->SetVarRange(2,0.,1.0001);// prob
 571   fQA->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
 572   fQA->SetVarRange(4,-0.5,1.5); // pid mask
 573   fQA->SetVarName(0,"centrality");
 574   fQA->SetVarName(1,"p_{t}");
 575   fQA->SetVarName(2,"prob");
 576   fQA->SetVarName(3,"#Psi");
 577   fQA->SetVarName(4,"PIDmask");
 578   if(fQAsw && fV2) fQA->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 579   if(fQAsw && fV2) fQA->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 580   if(fQAsw && fV2) fQA->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 581 //   if(fQAsw && fV2) fQA->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 582 //   fQA->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 583 //   fQA->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 584 //   fQA->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 585
 586   fQA2 = new AliFlowVZEROQA("v2CQA",5,binsQA);
 587   fQA2->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 588   fQA2->SetVarRange(1,0,7);  // pt
 589   fQA2->SetVarRange(2,0.,1.0001);// prob
 590   fQA2->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
 591   fQA2->SetVarRange(4,-0.5,1.5); // pid mask
 592   fQA2->SetVarName(0,"centrality");
 593   fQA2->SetVarName(1,"p_{t}");
 594   fQA2->SetVarName(2,"prob");
 595   fQA2->SetVarName(3,"#Psi");
 596   fQA2->SetVarName(4,"PIDmask");
 597   if(fQAsw && fV2) fQA2->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 598   if(fQAsw && fV2) fQA2->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 599   if(fQAsw && fV2) fQA2->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 600 //   if(fQAsw && fV2) fQA2->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 601 //   fQA2->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 602 //   fQA2->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 603 //   fQA2->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 604
 605   fQAv3 = new AliFlowVZEROQA("v3AQA",5,binsQAv3);
 606   fQAv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 607   fQAv3->SetVarRange(1,0,7);  // pt
 608   fQAv3->SetVarRange(2,0.,1.0001);// prob
 609   fQAv3->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
 610   fQAv3->SetVarRange(4,-0.5,1.5); // pid mask
 611   fQAv3->SetVarName(0,"centrality");
 612   fQAv3->SetVarName(1,"p_{t}");
 613   fQAv3->SetVarName(2,"prob");
 614   fQAv3->SetVarName(3,"#Psi");
 615   fQAv3->SetVarName(4,"PIDmask");
 616   if(fQAsw && fV3) fQAv3->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 617 //   if(fQAsw && fV3) fQAv3->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 618 //   if(fQAsw && fV3) fQAv3->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 619 //   if(fQAsw && fV2) fQAv3->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 620 //   fQAv3->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 621 //   fQAv3->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 622 //   fQAv3->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 623
 624   fQA2v3 = new AliFlowVZEROQA("v3CQA",5,binsQAv3);
 625   fQA2v3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 626   fQA2v3->SetVarRange(1,0,7);  // pt
 627   fQA2v3->SetVarRange(2,0.,1.0001);// prob
 628   fQA2v3->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
 629   fQA2v3->SetVarRange(4,-0.5,1.5); // pid mask
 630   fQA2v3->SetVarName(0,"centrality");
 631   fQA2v3->SetVarName(1,"p_{t}");
 632   fQA2v3->SetVarName(2,"prob");
 633   fQA2v3->SetVarName(3,"#Psi");
 634   fQA2v3->SetVarName(4,"PIDmask");
 635   if(fQAsw && fV3) fQA2v3->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 636 //   if(fQAsw && fV3) fQA2v3->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 637 //   if(fQAsw && fV3) fQA2v3->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 638 //   if(fQAsw && fV2) fQA2v3->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 639 //   fQA2v3->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 640 //   fQA2v3->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 641 //   fQA2v3->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 642
 643   fList->Add(fPhiRPv0A);
 644   fList->Add(fPhiRPv0C);
 645
 646   if(fQAsw && fV2){
 647     fList4->Add(fQA);
 648     fList4->Add(fQA2);
 649   }
 650
 651   fList2->Add(fPhiRPv0Av3);
 652   fList2->Add(fPhiRPv0Cv3);
 653
 654   if(fQAsw && fV3){
 655    fList4->Add(fQAv3);
 656    fList4->Add(fQA2v3);
 657   }
 658
 659   //  fList->Add(fTree); // comment if not needed
 660
 661   const Int_t nDCA = 300;
 662   Double_t DCAbin[nDCA+1];
 663   for(Int_t i=0;i <= nDCA;i++){
 664     DCAbin[i] = -3 +i*6.0/nDCA;
 665   }
 666
 667   char nameHistos[100];
 668   for(Int_t iC=0;iC < nCentrBin;iC++){
 669     snprintf(nameHistos,100,"fHdcaPtPiCent%i",iC);
 670     fHdcaPt[iC][0] = new TH2D(nameHistos,"DCA_{xy} for #pi;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 671     snprintf(nameHistos,100,"fHdcaPtKaCent%i",iC);
 672     fHdcaPt[iC][1] = new TH2D(nameHistos,"DCA_{xy} for K;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 673     snprintf(nameHistos,100,"fHdcaPtPrCent%i",iC);
 674     fHdcaPt[iC][2] = new TH2D(nameHistos,"DCA_{xy} for #bar{p};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 675     snprintf(nameHistos,100,"fHdcaPtElCent%i",iC);
 676     fHdcaPt[iC][3] = new TH2D(nameHistos,"DCA_{xy} for e;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 677     snprintf(nameHistos,100,"fHdcaPtDeCent%i",iC);
 678     fHdcaPt[iC][4] = new TH2D(nameHistos,"DCA_{xy} for #bar{d};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 679     snprintf(nameHistos,100,"fHdcaPtTrCent%i",iC);
 680     fHdcaPt[iC][5] = new TH2D(nameHistos,"DCA_{xy} for #bar{t};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 681     snprintf(nameHistos,100,"fHdcaPtHeCent%i",iC);
 682     fHdcaPt[iC][6] = new TH2D(nameHistos,"DCA_{xy} for #bar{^{3}He};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 683   }
 684
 685   if(fFillDCA && fQAsw){
 686     for(Int_t i=0;i<7;i++)
 687       for(Int_t iC=0;iC < nCentrBin;iC++)
 688         fList4->Add(fHdcaPt[iC][i]);
 689   }
 690   if(fIsMC){
 691     for(Int_t iC=0;iC < nCentrBin;iC++){
 692       snprintf(nameHistos,100,"fHdcaPtPiSecCent%i",iC);
 693       fHdcaPtSec[iC][0] = new TH2D(nameHistos,"DCA_{xy} for secondary #pi;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 694       snprintf(nameHistos,100,"fHdcaPtKaSecCent%i",iC);
 695       fHdcaPtSec[iC][1] = new TH2D(nameHistos,"DCA_{xy} for secondary K;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 696       snprintf(nameHistos,100,"fHdcaPtPrSecCent%i",iC);
 697       fHdcaPtSec[iC][2] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{p};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 698       snprintf(nameHistos,100,"fHdcaPtElSecCent%i",iC);
 699       fHdcaPtSec[iC][3] = new TH2D(nameHistos,"DCA_{xy} for secondary e;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 700       snprintf(nameHistos,100,"fHdcaPtDeSecCent%i",iC);
 701       fHdcaPtSec[iC][4] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{d};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 702       snprintf(nameHistos,100,"fHdcaPtTrSecCent%i",iC);
 703       fHdcaPtSec[iC][5] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{t};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 704       snprintf(nameHistos,100,"fHdcaPtHeSecCent%i",iC);
 705       fHdcaPtSec[iC][6] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{^{3}He};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
 706     }
 707
 708     if(fFillDCA && fQAsw){
 709       for(Int_t i=0;i<7;i++)
 710         for(Int_t iC=0;iC < nCentrBin;iC++)
 711           fList4->Add(fHdcaPtSec[iC][i]);
 712     }
 713   }
 714
 715   // Add TProfile Extra QA
 716   const Int_t nBinQApid = 2;
 717   Int_t binQApid[nBinQApid] = {nCentrTOF,200};
 718   const Int_t nbinsigma = 100;
 719   Double_t nsigmaQA[nbinsigma+1];
 720   for(Int_t i=0;i<nbinsigma+1;i++){
 721     nsigmaQA[i] = -10 + 20.0*i/nbinsigma;
 722   }
 723   fContQApid = new AliFlowVZEROResults("qaPID",nBinQApid,binQApid);
 724   fContQApid->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
 725   fContQApid->SetVarRange(1,0,20);  // charge
 726   fContQApid->SetVarName(0,"centrality");
 727   fContQApid->SetVarName(1,"p_{t}");
 728   fContQApid->AddSpecies("piTPC",nbinsigma,nsigmaQA);
 729   fContQApid->AddSpecies("piTOF",nbinsigma,nsigmaQA);
 730   fContQApid->AddSpecies("kaTPC",nbinsigma,nsigmaQA);
 731   fContQApid->AddSpecies("kaTOF",nbinsigma,nsigmaQA);
 732   fContQApid->AddSpecies("prTPC",nbinsigma,nsigmaQA);
 733   fContQApid->AddSpecies("prTOF",nbinsigma,nsigmaQA);
 734   if(fV2) fList->Add(fContQApid);
 735
 736   printf("Output creation ok!!\n\n\n\n");
 737
 738   fList->Add(fHctauPtEP);
 739   fList->Add(fHctauAt1EP);
 740
 741   fHKsPhi = new TH2D("hKsPhi","K^{0}_{s} #phi distributuion;v_{z} (cm);#phi (rad)",20,-10,10,20,0,2*TMath::Pi());
 742   fList->Add(fHKsPhi);
 743   fHKsPhiEP = new TH2D("hKsPhiEP","EP V0C #phi distributuion;v_{z} (cm);#phi (rad)",20,-10,10,20,0,2*TMath::Pi());
 744   fList->Add(fHKsPhiEP);
 745
 746   fHK0sMass = new TH2D("hK0sMass","K^{0}_{s} mass;p_{T} (GeV/#it{c});mass (GeV/#it{c}^{2})",20,0,5,400,0,1);
 747   fList->Add(fHK0sMass);
 748   fHK0sMass2 = new TH2D("hK0sMass2","K^{0}_{s} mass using secondary vertex;p_{T} (GeV/#it{c});mass (GeV/#it{c}^{2})",20,0,5,400,0,1);
 749   fList->Add(fHK0sMass2);
 750
 751   fHK0vsLambda= new TH2D("hK0vsLambda",";K^{0} mass;#Lambda mass",100,0,1,100,0.5,1.5);
 752   fList->Add(fHK0vsLambda);
 753
 754   // Post output data.
 755   if(fV2) PostData(1, fList);
 756   if(fV3) PostData(2, fList2);
 757   if(fIsMC) PostData(3, fList3);
 758   if(fQAsw) PostData(4, fList4);
 759
 760 }
 761
 762 //______________________________________________________________________________
 763 void AliAnalysisTaskVnV0::UserExec(Option_t *)
 764 {
 765     // Main loop
 766     // Called for each event
 767
 768     fgIsPsiComputed = kFALSE;
 769     fgPsi2v0a=999.;
 770     fgPsi2v0c=999.;
 771     fgPsi2tpc=999.;
 772     fgPsi3v0a=999.;
 773     fgPsi3v0c=999.;
 774     fgPsi3tpc=999.;
 775     fgPsi2v0aMC=999.;
 776     fgPsi2v0cMC=999.;
 777     fgPsi2tpcMC=999.;
 778     fgPsi3v0aMC=999.;
 779     fgPsi3v0cMC=999.;
 780     fgPsi3tpcMC=999.;
 781
 782     fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 783     if(!fOutputAOD){
 784         Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 785         this->Dump();
 786         return;
 787     }
 788
 789     Int_t run = fOutputAOD->GetRunNumber();
 790
 791     if(run != fRun){
 792         // Load the calibrations run dependent
 793       if(! fIsAfter2011) OpenInfoCalbration(run);
 794       fRun=run;
 795     }
 796
 797     fZvtx = GetVertex(fOutputAOD);
 798
 799
 800
 801     //Get the MC object
 802     if(fIsMC){
 803       AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 804       if (!mcHeader) {
 805         AliError("Could not find MC Header in AOD");
 806         return;
 807       }
 808     }
 809
 810     /*
 811     AliMCEvent* mcEvent = MCEvent();
 812     if (!mcEvent) {
 813       Printf("ERROR: Could not retrieve MC event");
 814       return;
 815     }
 816
 817     Double_t gReactionPlane = -999., gImpactParameter = -999.;
 818     //Get the MC header
 819     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
 820     if (headerH) {
 821       //Printf("=====================================================");
 822       //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());
 823       //Printf("=====================================================");
 824       gReactionPlane = headerH->ReactionPlaneAngle();
 825       gImpactParameter = headerH->ImpactParameter();
 826     }
 827
 828 */
 829
 830     if (TMath::Abs(fZvtx) < fVtxCut) {
 831       //Centrality
 832       Float_t v0Centr  = -10.;
 833       Float_t trkCentr  = -10.;
 834       AliCentrality *centrality = fOutputAOD->GetCentrality();
 835       if (centrality){
 836 //      printf("v0centr = %f -- tpccnetr%f\n",centrality->GetCentralityPercentile("V0M"),centrality->GetCentralityPercentile("TRK"));
 837         v0Centr  = centrality->GetCentralityPercentile("V0M");
 838         trkCentr = centrality->GetCentralityPercentile("TRK");
 839         //changed
 840       }
 841
 842       if(TMath::Abs(v0Centr - trkCentr) < 5.0 && v0Centr > 0){ // consistency cut on centrality selection
 843         fPID->SetDetResponse(fOutputAOD, v0Centr); // Set the PID object for each event!!!!
 844         Analyze(fOutputAOD,v0Centr); // Do analysis!!!!
 845
 846         fCentrality = v0Centr;
 847         if(fV2) fTree->Fill();
 848       }
 849     }
 850
 851 }
 852
 853 //________________________________________________________________________
 854 void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
 855 {
 856
 857   Int_t nusedForK0s=0;
 858   AliAODTrack *usedForK0s1[1000];
 859   AliAODTrack *usedForK0s2[1000];
 860
 861   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};
 862
 863   // Event plane resolution for v2
 864   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
 865                          0.446480,0.612705,0.712222,0.736200,0.697907,0.610114,0.481009,0.327402,0.182277};// V0C vs. centrality
 866
 867   Int_t iC = -1;
 868   if (v0Centr < 80){ // analysis only for 0-80% centrality classes
 869     //  if (v0Centr >0 && v0Centr < 80){ // analysis only for 0-80% centrality classes
 870     // changed
 871     fgIsPsiComputed = kTRUE;
 872
 873     // centrality bins
 874     if(v0Centr < 5) iC = 0;
 875     else if(v0Centr < 10) iC = 1;
 876     else if(v0Centr < 20) iC = 2;
 877     else if(v0Centr < 30) iC = 3;
 878     else if(v0Centr < 40) iC = 4;
 879     else if(v0Centr < 50) iC = 5;
 880     else if(v0Centr < 60) iC = 6;
 881     else if(v0Centr < 70) iC = 7;
 882     else iC = 8;
 883
 884     Int_t iCcal = iC;
 885
 886 /*
 887     if(nCentrBin==16){
 888        iC = v0Centr/5;
 889        if(iC >= nCentrBin) iC = nCentrBin-1;
 890     }
 891
 892     // centrality bins
 893     // changed
 894     if(v0Centr < 10 + 10./9) iC = 0;
 895     else if(v0Centr < 10 + 20./9) iC = 1;
 896     else if(v0Centr < 10 + 30./9) iC = 2;
 897     else if(v0Centr < 10 + 40./9) iC = 3;
 898     else if(v0Centr < 10 + 50./9) iC = 4;
 899     else if(v0Centr < 10 + 60./9) iC = 5;
 900     else if(v0Centr < 10 + 70./9) iC = 6;
 901     else if(v0Centr < 10 + 90./9) iC = 7;
 902     else if(v0Centr < 10 + 100./9) iC = 8;
 903     else if(v0Centr < 10 + 110./9) iC = 9;
 904     else if(v0Centr < 10 + 120./9) iC = 10;
 905     else if(v0Centr < 10 + 130./9) iC = 11;
 906     else if(v0Centr < 10 + 140./9) iC = 12;
 907     else if(v0Centr < 10 + 150./9) iC = 13;
 908     else if(v0Centr < 10 + 160./9) iC = 14;
 909     else if(v0Centr < 10 + 170./9) iC = 15;
 910     else iC = 16;
 911     if(iC >= nCentrBin) iC= nCentrBin - 1;
 912 */
 913
 914     //reset Q vector info
 915     Double_t Qxa2 = 0, Qya2 = 0;
 916     Double_t Qxc2 = 0, Qyc2 = 0;
 917     Double_t Qxa3 = 0, Qya3 = 0;
 918     Double_t Qxc3 = 0, Qyc3 = 0;
 919
 920     Int_t nAODTracks = aodEvent->GetNumberOfTracks();
 921
 922     AliAODMCHeader *mcHeader = NULL;
 923     TClonesArray *mcArray = NULL;
 924     Float_t evplaneMC = 0;
 925     if(fIsMC){
 926       mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 927
 928       if (mcHeader) {
 929         evplaneMC = mcHeader->GetReactionPlaneAngle();
 930         if(evplaneMC > TMath::Pi()/2 && evplaneMC <=  TMath::Pi()*3/2) evplaneMC-=TMath::Pi();
 931         else if(evplaneMC > TMath::Pi()*3/2) evplaneMC-=2*TMath::Pi();
 932         mcArray = (TClonesArray*)fOutputAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 933
 934         if(mcArray){
 935           Float_t QxMCv2[3] = {0,0,0};
 936           Float_t QyMCv2[3] = {0,0,0};
 937           Float_t QxMCv3[3] = {0,0,0};
 938           Float_t QyMCv3[3] = {0,0,0};
 939           Float_t EvPlaneMCV2[3] = {0,0,0};
 940           Float_t EvPlaneMCV3[3] = {0,0,0};
 941           Float_t etaMin[3] = {2.8,-3.6,-0.8}; // A-side, C-side M-barrel
 942           Float_t etaMax[3] = {4.88,-1.8,0.8};
 943
 944           // analysis on MC tracks
 945           Int_t nMCtrack = mcArray->GetEntries() ;
 946
 947           // EP computation with MC tracks
 948           for(Int_t iT=0;iT < nMCtrack;iT++){
 949             AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
 950             if(!mctr || !(mctr->IsPrimary()) || !(mctr->Charge()) || mctr->Pt() < 0.2) continue;
 951
 952             Float_t eta = mctr->Eta();
 953
 954             for(Int_t iD=0;iD<3;iD++){
 955               if(eta > etaMin[iD] && eta < etaMax[iD]){
 956                 Float_t phi = mctr->Phi();
 957                 QxMCv2[iD] += TMath::Cos(2*phi);
 958                 QyMCv2[iD] += TMath::Sin(2*phi);
 959                 QxMCv3[iD] += TMath::Cos(3*phi);
 960                 QyMCv3[iD] += TMath::Sin(3*phi);
 961               }
 962             }
 963           }
 964
 965           if(fV2){
 966             EvPlaneMCV2[0] = TMath::ATan2(QyMCv2[0],QxMCv2[0])/2.;
 967             EvPlaneMCV2[1] = TMath::ATan2(QyMCv2[1],QxMCv2[1])/2.;
 968             EvPlaneMCV2[2] = TMath::ATan2(QyMCv2[2],QxMCv2[2])/2.;
 969             fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
 970             fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
 971             fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
 972             fgPsi2v0aMC = EvPlaneMCV2[0];
 973             fgPsi2v0cMC = EvPlaneMCV2[1];
 974             fgPsi2tpcMC = EvPlaneMCV2[2];
 975           }
 976           if(fV3){
 977             EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
 978             EvPlaneMCV3[1] = TMath::ATan2(QyMCv3[1],QxMCv3[1])/3.;
 979             EvPlaneMCV3[2] = TMath::ATan2(QyMCv3[2],QxMCv3[2])/3.;
 980             fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
 981             fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
 982             fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
 983             fgPsi3v0aMC = EvPlaneMCV3[0];
 984             fgPsi3v0cMC = EvPlaneMCV3[1];
 985             fgPsi3tpcMC = EvPlaneMCV3[2];
 986           }
 987
 988           // flow A and C side
 989           Float_t xMCepAv2[5] = {iC,0/*charge*/,1,EvPlaneMCV2[0],1};
 990           Float_t xMCepCv2[5] = {iC,0/*charge*/,1,EvPlaneMCV2[1],1};
 991           Float_t xMCepAv3[5] = {iC,0/*charge*/,1,EvPlaneMCV3[0],1};
 992           Float_t xMCepCv3[5] = {iC,0/*charge*/,1,EvPlaneMCV3[1],1};
 993
 994           for(Int_t iT=0;iT < nMCtrack;iT++){
 995             AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
 996             if(!mctr || !(mctr->IsPhysicalPrimary()) || !(mctr->Charge()) || TMath::Abs(mctr->Eta()) > 0.8 || mctr->Pt() < 0.2) continue;
 997             Int_t iS = TMath::Abs(mctr->GetPdgCode());
 998             Int_t charge = mctr->Charge();
 999             Float_t pt = mctr->Pt();
1000             Float_t phi = mctr->Phi();
1001
1002             if(charge > 0){
1003               xMCepAv2[1] = 1;
1004               xMCepCv2[1] = 1;
1005               xMCepAv3[1] = 1;
1006               xMCepCv3[1] = 1;
1007             }
1008             else{
1009               xMCepAv2[1] = -1;
1010               xMCepCv2[1] = -1;
1011               xMCepAv3[1] = -1;
1012               xMCepCv3[1] = -1;
1013             }
1014
1015             fContAllChargesMCA->Fill(0,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1016             fContAllChargesMCC->Fill(0,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1017             fContAllChargesMCAv3->Fill(0,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1018             fContAllChargesMCCv3->Fill(0,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
1019
1020             if(iS==11){
1021               fContAllChargesMCA->Fill(4,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1022               fContAllChargesMCC->Fill(4,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1023               fContAllChargesMCAv3->Fill(4,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1024               fContAllChargesMCCv3->Fill(4,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
1025             }
1026             else if(iS==13){
1027               fContAllChargesMCA->Fill(5,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1028               fContAllChargesMCC->Fill(5,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1029               fContAllChargesMCAv3->Fill(5,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1030               fContAllChargesMCCv3->Fill(5,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
1031             }
1032             else if(iS==211){
1033               fContAllChargesMCA->Fill(1,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1034               fContAllChargesMCC->Fill(1,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1035               fContAllChargesMCAv3->Fill(1,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1036               fContAllChargesMCCv3->Fill(1,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
1037             }
1038             else if(iS==321){
1039               fContAllChargesMCA->Fill(2,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1040               fContAllChargesMCC->Fill(2,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1041               fContAllChargesMCAv3->Fill(2,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1042               fContAllChargesMCCv3->Fill(2,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
1043             }
1044             else if(iS==2212){
1045               fContAllChargesMCA->Fill(3,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1046               fContAllChargesMCC->Fill(3,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1047               fContAllChargesMCAv3->Fill(3,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1048               fContAllChargesMCCv3->Fill(3,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
1049             }
1050           }
1051         }
1052       }
1053     }
1054
1055     // TPC EP needed for resolution studies (TPC subevent)
1056     Double_t Qx2 = 0, Qy2 = 0;
1057     Double_t Qx3 = 0, Qy3 = 0;
1058
1059     for(Int_t iT = 0; iT < nAODTracks; iT++) {
1060
1061       AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
1062
1063       if (!aodTrack){
1064         aodTrack->Delete();
1065         continue;
1066       }
1067
1068       Bool_t trkFlag = aodTrack->TestFilterBit(1);
1069
1070       if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag)
1071         continue;
1072
1073       Double_t b[2] = {-99., -99.};
1074       Double_t bCov[3] = {-99., -99., -99.};
1075       if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
1076         continue;
1077
1078       if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
1079         continue;
1080
1081       Qx2 += TMath::Cos(2*aodTrack->Phi());
1082       Qy2 += TMath::Sin(2*aodTrack->Phi());
1083       Qx3 += TMath::Cos(3*aodTrack->Phi());
1084       Qy3 += TMath::Sin(3*aodTrack->Phi());
1085
1086     }
1087
1088     evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
1089     evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
1090
1091     fgPsi2tpc = evPlAng2;
1092     fgPsi3tpc = evPlAng3;
1093
1094     SelectK0s();
1095
1096     //V0 info
1097     AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
1098
1099     for (Int_t iv0 = 0; iv0 < 64; iv0++) {
1100       Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
1101       Float_t multv0 = aodV0->GetMultiplicity(iv0);
1102
1103       if(! fIsAfter2011){
1104         if(! fIsMC){
1105           if (iv0 < 32){ // V0C
1106             Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1107             Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1108             Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1109             Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1110           } else {       // V0A
1111             Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1112             Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1113             Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1114             Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1115           }
1116         }
1117         else{
1118           if (iv0 < 32){ // V0C
1119             Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1120             Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1121             Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1122             Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1123           } else {       // V0A
1124             Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1125             Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1126             Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1127             Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1128           }
1129         }
1130       }
1131     }
1132
1133     //grab for each centrality the proper histo with the Qx and Qy to do the recentering
1134     Double_t Qxamean2 = fMeanQ[iCcal][1][0];
1135     Double_t Qxarms2  = fWidthQ[iCcal][1][0];
1136     Double_t Qyamean2 = fMeanQ[iCcal][1][1];
1137     Double_t Qyarms2  = fWidthQ[iCcal][1][1];
1138     Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
1139     Double_t Qxarms3  = fWidthQv3[iCcal][1][0];
1140     Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
1141     Double_t Qyarms3  = fWidthQv3[iCcal][1][1];
1142
1143     Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
1144     Double_t Qxcrms2  = fWidthQ[iCcal][0][0];
1145     Double_t Qycmean2 = fMeanQ[iCcal][0][1];
1146     Double_t Qycrms2  = fWidthQ[iCcal][0][1];
1147     Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
1148     Double_t Qxcrms3  = fWidthQv3[iCcal][0][0];
1149     Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
1150     Double_t Qycrms3  = fWidthQv3[iCcal][0][1];
1151
1152     Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
1153     Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
1154     Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
1155     Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
1156     Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
1157     Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
1158     Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
1159     Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
1160
1161     if(! fIsAfter2011){
1162       if(! fIsMC){
1163         evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
1164         evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
1165         evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
1166         evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
1167       }
1168       else{
1169         evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
1170         evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
1171         evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
1172         evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
1173       }
1174     }
1175     else{
1176       AliEventplane *ep =  aodEvent->GetEventplane();
1177       evPlAngV0ACor2 = ep->GetEventplane("V0A", aodEvent, 2);
1178       evPlAngV0CCor2 = ep->GetEventplane("V0C", aodEvent, 2);
1179       evPlAngV0ACor3 = ep->GetEventplane("V0A", aodEvent, 3);
1180       evPlAngV0CCor3 = ep->GetEventplane("V0C", aodEvent, 3);
1181     }
1182
1183
1184     fgPsi2v0a = evPlAngV0ACor2;
1185     fgPsi2v0c = evPlAngV0CCor2;
1186     fgPsi3v0a = evPlAngV0ACor3;
1187     fgPsi3v0c = evPlAngV0CCor3;
1188
1189     fHKsPhiEP->Fill(fZvtx,fgPsi2v0c);
1190     //loop track and get pid
1191     for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
1192       AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
1193
1194       if (!aodTrack){
1195         aodTrack->Delete();
1196         continue;
1197       }
1198
1199       Bool_t trkFlag = aodTrack->TestFilterBit(1); // TPC only tracks
1200       if(fFillDCA)
1201         trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
1202
1203       if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
1204         continue;
1205       }
1206
1207       Double_t b[2] = {-99., -99.};
1208       Double_t bCov[3] = {-99., -99., -99.};
1209       if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
1210         continue;
1211
1212       if (!fFillDCA && ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4)))
1213         continue;
1214
1215       if(fFillDCA && (TMath::Abs(b[0]) > 3.0 || TMath::Abs(b[1]) > 3))
1216         continue;
1217
1218       // re-map the container in an array to do the analysis for V0A and V0C within a loop
1219       Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
1220       AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
1221       AliFlowVZEROQA *QA[2] = {fQA,fQA2};
1222
1223       Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
1224       AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
1225       AliFlowVZEROQA *QAv3[2] = {fQAv3,fQA2v3};
1226
1227       // Fill MC results
1228       if(fIsMC && mcArray){
1229         fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1230         Float_t tofMismProbMC = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1231
1232         Float_t xMC[5] = {iC,aodTrack->Charge(),1,evplaneMC,fPID->GetCurrentMask(1)&&tofMismProbMC < 0.5}; // to fill analysis v2 container
1233
1234         Float_t v2mc = TMath::Cos(2*(aodTrack->Phi() - evplaneMC));
1235
1236         fContAllChargesMC->Fill(0,aodTrack->Pt(),v2mc,xMC);
1237
1238         Int_t iS = TMath::Abs(((AliAODMCParticle*)mcArray->At(TMath::Abs(aodTrack->GetLabel())))->GetPdgCode());
1239         if(iS==11){
1240           fContAllChargesMC->Fill(4,aodTrack->Pt(),v2mc,xMC);
1241         }
1242         else if(iS==13){
1243           fContAllChargesMC->Fill(5,aodTrack->Pt(),v2mc,xMC);
1244         }
1245         else if(iS==211){
1246           fContAllChargesMC->Fill(1,aodTrack->Pt(),v2mc,xMC);
1247         }
1248         else if(iS==321){
1249           fContAllChargesMC->Fill(2,aodTrack->Pt(),v2mc,xMC);
1250         }
1251         else if(iS==2212){
1252           fContAllChargesMC->Fill(3,aodTrack->Pt(),v2mc,xMC);
1253         }
1254       }
1255
1256       for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1257
1258         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)
1259
1260         Float_t v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1261         Float_t v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1262
1263         fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1264         Float_t dedx = fPID->GetDeDx();//aodTrack->GetTPCsignal();
1265         Float_t *probRead = fPID->GetProb();
1266         Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1267         Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1268         Float_t x[6] = {iC,aodTrack->Charge(),1,evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v2 container
1269         Float_t x3[6] = {iC,aodTrack->Charge(),1,evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v3 container
1270
1271         // in case fill DCA info
1272         if(fFillDCA){
1273           if(TMath::Abs(b[0]) > 0.1){
1274             x[5] = 1;
1275             x3[5] = 1;
1276           }
1277           if(TMath::Abs(b[0]) > 0.3){
1278             x[5] = 2;
1279             x3[5] = 2;
1280           }
1281           if(fIsMC && mcArray){
1282             if(!((AliAODMCParticle*)mcArray->At(TMath::Abs(aodTrack->GetLabel())))->IsPhysicalPrimary()){
1283               x[5] += 3;
1284               x3[5] += 3;
1285             }
1286           }
1287         }
1288
1289         // Fill no PID
1290         if(fV2) contV0[iV0]->Fill(0,aodTrack->Pt(),v2V0,x);
1291         if(fV3) contV0v3[iV0]->Fill(0,aodTrack->Pt(),v3V0,x3);
1292
1293
1294         Double_t dedxExp[8];
1295         Float_t tof = -1;
1296         Double_t inttimes[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
1297         Double_t expTOFsigma[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
1298
1299         Float_t nsigmaTPC[8];
1300         Float_t nsigmaTOF[8];
1301
1302         if(aodTrack->GetDetPid()){ // check the PID object is available
1303           for(Int_t iS=0;iS < 8;iS++){
1304             dedxExp[iS] = fPID->GetExpDeDx(aodTrack,iS);
1305             nsigmaTPC[iS] = (dedx - dedxExp[iS])/(dedxExp[iS]*0.07);
1306             //      printf("TPC %i = %f (%f %f)\n",iS, nsigmaTPC[iS],dedx, dedxExp[iS]);
1307           }
1308
1309           if(fPID->GetCurrentMask(1)){ // if TOF is present
1310             Float_t ptrack = aodTrack->P();
1311             tof = aodTrack->GetTOFsignal() - fPID->GetESDpid()->GetTOFResponse().GetStartTime(ptrack);
1312             aodTrack->GetIntegratedTimes(inttimes);
1313
1314             for(Int_t iS=5;iS < 8;iS++) // extra info for light nuclei
1315               inttimes[iS] = inttimes[0] / ptrack * mass[iS] * TMath::Sqrt(1+ptrack*ptrack/mass[iS]/mass[iS]);
1316
1317             for(Int_t iS=0;iS<8;iS++){
1318               expTOFsigma[iS] = fPID->GetESDpid()->GetTOFResponse().GetExpectedSigma(ptrack, inttimes[iS], mass[iS]);
1319               nsigmaTOF[iS] = (tof - inttimes[iS])/expTOFsigma[iS];
1320               //              printf("TOF %i = %f\n",iS, nsigmaTOF[iS]);
1321             }
1322           }
1323         }
1324
1325         Float_t deltaPhiV0 = aodTrack->Phi() - evPlAngV0[iV0];
1326         if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi();
1327         else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi();
1328         if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi();
1329         else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi();
1330
1331         Float_t deltaPhiV0v3 = aodTrack->Phi() - evPlAngV0v3[iV0];
1332         if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi();
1333         else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi();
1334         if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi();
1335         else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi();
1336
1337         // variable to fill QA container
1338         Float_t xQA[5] = {iC,aodTrack->Pt(), 0.0,deltaPhiV0,x[4]}; // v2
1339         Float_t xQA3[5] = {iC,aodTrack->Pt(), 0.0,deltaPhiV0v3,x[4]}; // v3
1340
1341         // extra QA TProfiles
1342         if(iV0==1 && aodTrack->Pt() < 20 && fPID->GetCurrentMask(0) && fPID->GetCurrentMask(1)){
1343           Float_t xQApid[2] = {iC,aodTrack->Pt()};
1344           fContQApid->Fill(0,nsigmaTPC[2],v2V0,xQApid); // v2 TPC (V0C) w.r.t pions
1345           fContQApid->Fill(1,nsigmaTOF[2],v2V0,xQApid); // v2 TOF (V0C) w.r.t. pions
1346           fContQApid->Fill(2,nsigmaTPC[3],v2V0,xQApid); // v2 TPC (V0C) w.r.t kaons
1347           fContQApid->Fill(3,nsigmaTOF[3],v2V0,xQApid); // v2 TOF (V0C) w.r.t. kaons
1348           fContQApid->Fill(4,nsigmaTPC[4],v2V0,xQApid); // v2 TPC (V0C) w.r.t protons
1349           fContQApid->Fill(5,nsigmaTOF[4],v2V0,xQApid); // v2 TOF (V0C) w.r.t. protons
1350         }
1351
1352         // QA fill
1353         if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid() || dedx < 10. || aodTrack->Pt() < 0 || aodTrack->Pt() > 7){}
1354         else{
1355           if(TMath::Abs(nsigmaTPC[2])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[2])<5))){ //pi
1356             xQA[2] = prob[2];
1357             xQA3[2] = xQA[2];
1358             if(fQAsw && fV2) QA[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA);
1359             if(fQAsw && fV3) QAv3[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA3);
1360           }
1361           if(TMath::Abs(nsigmaTPC[3])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[3])<5))){ //K
1362             xQA[2] = prob[3];
1363             xQA3[2] = xQA[2];
1364             if(fQAsw && fV2) QA[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA);
1365 //          if(fQAsw && fV3) QAv3[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA3);
1366           }
1367           if(TMath::Abs(nsigmaTPC[4])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[4])<5))){//p
1368             xQA[2] = prob[4];
1369             xQA3[2] = xQA[2];
1370             if(fQAsw && fV2) QA[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA);
1371 //          if(fQAsw && fV3) QAv3[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA3);
1372           }
1373           if(TMath::Abs(nsigmaTPC[0])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[0])<5))){//e
1374             xQA[2] = prob[0];
1375             xQA3[2] = xQA[2];
1376 //          if(fQAsw && fV2) QA[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA);
1377 //          if(fQAsw && fV3) QAv3[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA3);
1378           }
1379           if(TMath::Abs(nsigmaTPC[5])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[5])<5))){//d
1380             xQA[2] = prob[5];
1381             xQA3[2] = xQA[2];
1382             //    if(fQAsw && fV2) QA[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA);
1383             //    if(fQAsw && fV3) QAv3[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA3);
1384           }
1385           if(TMath::Abs(nsigmaTPC[6])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[6])<5))){//t
1386             xQA[2] = prob[6];
1387             xQA3[2] = xQA[2];
1388             //    if(fQAsw && fV2) QA[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA);
1389             //    if(fQAsw && fV3) QAv3[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA3);
1390           }
1391           if(TMath::Abs(nsigmaTPC[7])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[7])<5))){//He3
1392             xQA[2] = prob[7];
1393             xQA3[2] = xQA[2];
1394             //    if(fQAsw && fV2) QA[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA);
1395             //    if(fQAsw && fV3) QAv3[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA3);
1396           }
1397         }
1398
1399         //pid selection
1400         if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID and PID object strictly required (very important!!!!)
1401         else if(prob[2] > 0.6){ // pi
1402           x[2] = prob[2];
1403           x3[2] = x[2];
1404           if(TMath::Abs(nsigmaTPC[2]) < 5){ // TPC 5 sigma extra cut to accept the track
1405             if(fV2) contV0[iV0]->Fill(1,aodTrack->Pt(),v2V0,x);
1406             if(fV3) contV0v3[iV0]->Fill(1,aodTrack->Pt(),v3V0,x3);
1407             if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][0]->Fill(aodTrack->Pt(),b[0]);
1408             else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][0]->Fill(aodTrack->Pt(),b[0]);
1409           }
1410         }
1411         else if(prob[3] > 0.6){ // K
1412           x[2] = prob[3];
1413           x3[2] = x[2];
1414           if(TMath::Abs(nsigmaTPC[3]) < 5){
1415             if(fV2) contV0[iV0]->Fill(2,aodTrack->Pt(),v2V0,x);
1416             if(fV3) contV0v3[iV0]->Fill(2,aodTrack->Pt(),v3V0,x3);
1417             if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][1]->Fill(aodTrack->Pt(),b[0]);
1418             else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][1]->Fill(aodTrack->Pt(),b[0]);
1419           }
1420         }
1421         else if(prob[4] > 0.6){ // p
1422           x[2] = prob[4];
1423           x3[2] = x[2];
1424           if(TMath::Abs(nsigmaTPC[4]) < 5){
1425             if(fV2) contV0[iV0]->Fill(3,aodTrack->Pt(),v2V0,x);
1426             if(fV3) contV0v3[iV0]->Fill(3,aodTrack->Pt(),v3V0,x3);
1427             if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][2]->Fill(aodTrack->Pt(),b[0]);
1428             else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][2]->Fill(aodTrack->Pt(),b[0]);
1429           }
1430         }
1431         else if(prob[0] > 0.6){ // e
1432           x[2] = prob[0];
1433           x3[2] = x[2];
1434           if(TMath::Abs(nsigmaTPC[0]) < 5){
1435             if(fV2) contV0[iV0]->Fill(4,aodTrack->Pt(),v2V0,x);
1436             if(fV3) contV0v3[iV0]->Fill(4,aodTrack->Pt(),v3V0,x3);
1437             if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][3]->Fill(aodTrack->Pt(),b[0]);
1438             else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][3]->Fill(aodTrack->Pt(),b[0]);
1439           }
1440         }
1441         else if(prob[1] > 0.6){ // mu
1442           x[2] = prob[1];
1443           x3[2] = x[2];
1444           if(TMath::Abs(nsigmaTPC[1]) < 5){
1445             if(fV2) contV0[iV0]->Fill(8,aodTrack->Pt(),v2V0,x);
1446             if(fV3) contV0v3[iV0]->Fill(8,aodTrack->Pt(),v3V0,x3);
1447           }
1448         }
1449         else if(prob[5] > 0.6){ // d
1450           x[2] = prob[5];
1451           x3[2] = x[2];
1452           if(TMath::Abs(nsigmaTPC[5]) < 5){
1453             if(fV2) contV0[iV0]->Fill(5,aodTrack->Pt(),v2V0,x);
1454             if(fV3) contV0v3[iV0]->Fill(5,aodTrack->Pt(),v3V0,x3);
1455             if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][4]->Fill(aodTrack->Pt(),b[0]);
1456             else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][4]->Fill(aodTrack->Pt(),b[0]);
1457           }
1458         }
1459         else if(prob[6] > 0.6){ // t
1460           x[2] = prob[6];
1461           x3[2] = x[2];
1462           if(TMath::Abs(nsigmaTPC[6]) < 5){
1463             if(fV2) contV0[iV0]->Fill(6,aodTrack->Pt(),v2V0,x);
1464             if(fV3) contV0v3[iV0]->Fill(6,aodTrack->Pt(),v3V0,x3);
1465             if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][5]->Fill(aodTrack->Pt(),b[0]);
1466             else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][5]->Fill(aodTrack->Pt(),b[0]);
1467           }
1468         }
1469         else if(prob[7] > 0.6){ // He3
1470           x[2] = prob[7];
1471           x3[2] = x[2];
1472           if(TMath::Abs(nsigmaTPC[7]) < 5){
1473             if(fV2) contV0[iV0]->Fill(7,aodTrack->Pt()*2,v2V0,x);
1474             if(fV3) contV0v3[iV0]->Fill(7,aodTrack->Pt()*2,v3V0,x3);
1475             if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][6]->Fill(aodTrack->Pt(),b[0]);
1476             else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][6]->Fill(aodTrack->Pt(),b[0]);
1477           }
1478         }
1479
1480         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)
1481           fPID->ResetDetOR(1); // exclude TOF from PID
1482           tofMismProb = 0;
1483
1484           fPID->ComputeProb(aodTrack,fOutputAOD);
1485           dedx = fPID->GetDeDx();//aodTrack->GetTPCsignal();
1486           probRead = fPID->GetProb();
1487
1488           fPID->SetDetOR(1); // include TOF for PID
1489         }
1490         Float_t probTPC[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]}; // TPC stand alone prbabilities
1491
1492         //pid selection TPC S.A. with TOF matching
1493         x[4]*=2; // set the mask to 2 id TOF is present
1494         if(x[4]<1 || !(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID S.A. PID in TOF acceptance
1495         else if(probTPC[2] > 0.6){ // pi
1496           x[2] = probTPC[2];
1497           x3[2] = x[2];
1498           if(TMath::Abs(nsigmaTPC[2]) < 5){
1499             if(fV2) contV0[iV0]->Fill(1,aodTrack->Pt(),v2V0,x);
1500             if(fV3) contV0v3[iV0]->Fill(1,aodTrack->Pt(),v3V0,x3);
1501           }
1502         }
1503         else if(probTPC[3] > 0.6){ // K
1504           x[2] = probTPC[3];
1505           x3[2] = x[2];
1506           if(TMath::Abs(nsigmaTPC[3]) < 5){
1507             if(fV2) contV0[iV0]->Fill(2,aodTrack->Pt(),v2V0,x);
1508             if(fV3) contV0v3[iV0]->Fill(2,aodTrack->Pt(),v3V0,x3);
1509           }
1510         }
1511         else if(probTPC[4] > 0.6){ // p
1512           x[2] = probTPC[4];
1513           x3[2] = x[2];
1514           if(TMath::Abs(nsigmaTPC[4]) < 5){
1515             if(fV2) contV0[iV0]->Fill(3,aodTrack->Pt(),v2V0,x);
1516             if(fV3) contV0v3[iV0]->Fill(3,aodTrack->Pt(),v3V0,x3);
1517           }
1518         }
1519         else if(probTPC[0] > 0.6){ // e
1520           x[2] = probTPC[0];
1521           x3[2] = x[2];
1522           if(TMath::Abs(nsigmaTPC[0]) < 5){
1523             if(fV2) contV0[iV0]->Fill(4,aodTrack->Pt(),v2V0,x);
1524             if(fV3) contV0v3[iV0]->Fill(4,aodTrack->Pt(),v3V0,x3);
1525           }
1526         }
1527         else if(probTPC[1] > 0.6){ // mu
1528           x[2] = probTPC[1];
1529           x3[2] = x[2];
1530           if(TMath::Abs(nsigmaTPC[1]) < 5){
1531             if(fV2) contV0[iV0]->Fill(8,aodTrack->Pt(),v2V0,x);
1532             if(fV3) contV0v3[iV0]->Fill(8,aodTrack->Pt(),v3V0,x3);
1533           }
1534         }
1535         else if(probTPC[5] > 0.6){ // d
1536           x[2] = probTPC[5];
1537           x3[2] = x[2];
1538           if(TMath::Abs(nsigmaTPC[5]) < 5){
1539             if(fV2) contV0[iV0]->Fill(5,aodTrack->Pt(),v2V0,x);
1540             if(fV3) contV0v3[iV0]->Fill(5,aodTrack->Pt(),v3V0,x3);
1541           }
1542         }
1543         else if(probTPC[6] > 0.6){ // t
1544           x[2] = probTPC[6];
1545           x3[2] = x[2];
1546           if(TMath::Abs(nsigmaTPC[6]) < 5){
1547             if(fV2) contV0[iV0]->Fill(6,aodTrack->Pt(),v2V0,x);
1548             if(fV3) contV0v3[iV0]->Fill(6,aodTrack->Pt(),v3V0,x3);
1549           }
1550         }
1551         else if(probTPC[7] > 0.6){ // He3
1552           x[2] = probTPC[7];
1553           x3[2] = x[2];
1554           if(TMath::Abs(nsigmaTPC[7]) < 5){
1555             if(fV2) contV0[iV0]->Fill(7,aodTrack->Pt()*2,v2V0,x);
1556             if(fV3) contV0v3[iV0]->Fill(7,aodTrack->Pt()*2,v3V0,x3);
1557           }
1558         }
1559       } // end side loop
1560     } // end track loop
1561
1562     // my V0 loop
1563     for(Int_t imy=0;imy<fNK0s;imy++){
1564       Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
1565       Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
1566
1567       AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
1568       AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
1569
1570       for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1571         Float_t x[6] = {iC,-1/*my K0s are negative for convention*/,1,evPlAngV0[iV0],1,0}; // to fill analysis v2 container
1572         Float_t x3[6] = {iC,-1,1,evPlAngV0v3[iV0],1,0}; // to fill analysis v3 container
1573
1574         Float_t v2V0 = TMath::Cos(2*(fPhiK0s[imy] - evPlAngV0[iV0]));
1575         Float_t v3V0 = TMath::Cos(3*(fPhiK0s[imy] - evPlAngV0v3[iV0]));
1576         if(fV2) contV0[iV0]->Fill(9,fPtK0s[imy],v2V0,x);
1577         if(fV3) contV0v3[iV0]->Fill(9,fPtK0s[imy],v3V0,x3);
1578       }
1579     }
1580
1581     // V0 loop
1582     Int_t nV0s = fOutputAOD->GetNumberOfV0s();
1583     AliAODv0 *myV0;
1584     Double_t dQT, dALPHA, dPT, dMASS=0.0;
1585     for (Int_t i=0; i!=nV0s; ++i) {
1586       myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
1587       if(!myV0) continue;
1588       if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > fEtaCut) continue; // skipping low momentum
1589       Int_t pass = PassesAODCuts(myV0,fOutputAOD,0);
1590       if(pass) {
1591         dMASS = myV0->MassK0Short();
1592         pass = 3;
1593         fHK0sMass2->Fill(myV0->Pt(),dMASS);
1594       }
1595       if(TMath::Abs(dMASS-0.497)/0.005 > 3){
1596         pass = PassesAODCuts(myV0,fOutputAOD,1);
1597         if(pass) dMASS = myV0->MassLambda();
1598         if(pass==2) dMASS = myV0->MassAntiLambda();
1599       }
1600
1601       if(pass){// 1 lambda, 2 antilambda, 3=K0s
1602         dPT=myV0->Pt();
1603         dQT=myV0->PtArmV0();
1604         dALPHA=myV0->AlphaV0();
1605
1606         Int_t iPos, iNeg;
1607         AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0);
1608         if(iT->Charge()>0) {
1609           iPos = 0; iNeg = 1;
1610         } else {
1611           iPos = 1; iNeg = 0;
1612         }
1613
1614         // check if one of the daugthers was already used
1615         if(pass == 3 && TMath::Abs(dMASS-0.497)/0.005 < 1){
1616           fHKsPhi->Fill(fZvtx, myV0->Phi());
1617         }
1618
1619         if(pass == 1000){ // disable
1620           Bool_t used = kFALSE;
1621           for(Int_t ii=0;ii<nusedForK0s;ii++){
1622             if(myV0->GetDaughter(iNeg) == usedForK0s1[ii] || myV0->GetDaughter(iPos) == usedForK0s2[ii]){
1623               used = kTRUE;
1624             }
1625           }
1626           if((!used) && nusedForK0s < 1000){
1627             nusedForK0s++;
1628             usedForK0s1[nusedForK0s] = (AliAODTrack *) myV0->GetDaughter(iNeg);
1629             usedForK0s2[nusedForK0s] = (AliAODTrack *) myV0->GetDaughter(iPos);
1630             printf("accepted\n");
1631           }
1632           else{
1633             dMASS = 0;
1634             printf("rejected\n");
1635           }
1636         }
1637
1638         iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1639         AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1640
1641         // re-map the container in an array to do the analysis for V0A and V0C within a loop
1642         Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
1643         AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
1644
1645         Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
1646         AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
1647
1648         for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1649
1650           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)
1651
1652           Float_t v2V0 = TMath::Cos(2*(myV0->Phi() - evPlAngV0[iV0]));
1653           Float_t v3V0 = TMath::Cos(3*(myV0->Phi() - evPlAngV0v3[iV0]));
1654
1655           Float_t deltaphi = myV0->Phi()- evPlAngV0[iV0];
1656           if(deltaphi > TMath::Pi()) deltaphi -= 2*TMath::Pi();
1657           if(deltaphi < -TMath::Pi()) deltaphi += 2*TMath::Pi();
1658
1659           Float_t x[6] = {iC,1,1,evPlAngV0[iV0],1,0}; // to fill analysis v2 container
1660           Float_t x3[6] = {iC,1,1,evPlAngV0v3[iV0],1,0}; // to fill analysis v3 container
1661
1662           Float_t decaylength = myV0->DecayLengthXY(fOutputAOD->GetPrimaryVertex());
1663           //      printf("decay length = %f\n",decaylength);
1664
1665           if(pass==2){ // anti-lambda charge = -1
1666             x[1] = -1;
1667             x3[1] = -1;
1668           }
1669
1670           if(decaylength < fMinDistV0) pass = 0;
1671           if(decaylength > fMaxDistV0) pass = 0;
1672
1673           Float_t nsigma = 0;
1674           if(pass < 3)
1675             nsigma = TMath::Abs(dMASS-1.116)/0.0016;
1676           else if(pass == 3)
1677             nsigma = TMath::Abs(dMASS-0.497)/0.005;
1678
1679           if(nsigma < 1)
1680             x[2] = 0.95;
1681           else if(nsigma < 2)
1682             x[2] = 0.85;
1683           else if(nsigma < 3)
1684             x[2] = 0.75;
1685           else if(nsigma < 4)
1686             x[2] = 0.65;
1687           else
1688             x[2] = 0.5;
1689
1690           x3[2] = x[2];
1691
1692           // Fill Container for lambda and Ks
1693           if(fV2 && pass == 3 && x[2] > 0.6){
1694             contV0[iV0]->Fill(9,myV0->Pt(),v2V0,x);
1695             fHctauPtEP->Fill(myV0->Pt(),deltaphi,decaylength);//ciao
1696             if(myV0->Pt() < 1.1 && myV0->Pt() > 0.9) fHctauAt1EP->Fill(decaylength,deltaphi);
1697           }
1698           if(fV3 && pass == 3 && x[2] > 0.6) contV0v3[iV0]->Fill(9,myV0->Pt(),v3V0,x3);
1699           if(fV2 && pass < 3 && x[2] > 0.6) contV0[iV0]->Fill(10,myV0->Pt(),v2V0,x);
1700           if(fV3 && pass < 3 && x[2] > 0.6) contV0v3[iV0]->Fill(10,myV0->Pt(),v3V0,x3);
1701
1702           if(pass < 3){ // lambda
1703             AliAODTrack* aodTrack = iT;
1704             if(pass==2) aodTrack=jT;
1705
1706             v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1707             v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1708
1709             fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1710             Float_t *probRead = fPID->GetProb();
1711             Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1712             Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1713
1714             if(prob[4] < 0.61) prob[4] = 0.61;
1715
1716             Float_t xdec[6] = {iC,aodTrack->Charge(),prob[4],evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v2 container
1717             Float_t xdec3[6] = {iC,aodTrack->Charge(),prob[4],evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v3 container
1718
1719             // Fill Container for (anti)proton from lambda
1720             if(nsigma < 2 && xdec[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
1721               if(fV2) contV0[iV0]->Fill(11,aodTrack->Pt(),v2V0,xdec);
1722               if(fV3) contV0v3[iV0]->Fill(11,aodTrack->Pt(),v3V0,xdec3);
1723             }
1724           }
1725           else if(pass == 3){
1726             AliAODTrack* aodTrack = iT;
1727
1728             v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1729             v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1730
1731             fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1732             Float_t *probRead = fPID->GetProb();
1733             Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1734             Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1735
1736             if(prob[2] < 0.61) prob[2] = 0.61;
1737
1738             Float_t xdec[6] = {iC,aodTrack->Charge(),prob[2],evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to
1739             Float_t xdec3[6] = {iC,aodTrack->Charge(),prob[2],evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to
1740
1741             if(nsigma < 2 && xdec[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
1742               if(fV2) contV0[iV0]->Fill(12,aodTrack->Pt(),v2V0,xdec);
1743               if(fV3) contV0v3[iV0]->Fill(12,aodTrack->Pt(),v3V0,xdec3);
1744             }
1745
1746             aodTrack = jT;
1747             v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1748             v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1749
1750             fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1751             Float_t *probRead2 = fPID->GetProb();
1752             Float_t prob2[8] = {probRead2[0],probRead2[1],probRead2[2],probRead2[3],probRead2[4],probRead2[5],probRead2[6],probRead2[7]};
1753             Float_t tofMismProb2 = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1754
1755             if(prob2[2] < 0.61) prob2[2] = 0.61;
1756
1757             Float_t xdecB[6] = {iC,aodTrack->Charge(),prob2[2],evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb2 < 0.5,0}; // to
1758             Float_t xdecB3[6] = {iC,aodTrack->Charge(),prob2[2],evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb2 < 0.5,0}; // to
1759
1760             if(nsigma < 2 && xdecB[2] > 0.6  && TMath::Abs(aodTrack->Eta()) < 0.8){
1761               if(fV2) contV0[iV0]->Fill(12,aodTrack->Pt(),v2V0,xdecB);
1762               if(fV3) contV0v3[iV0]->Fill(12,aodTrack->Pt(),v3V0,xdecB3);
1763             }
1764           }
1765         }
1766
1767       }
1768     } // end loop on V0
1769
1770
1771     // Fill EP distribution histograms
1772     if(fV2) fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
1773     if(fV2) fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
1774
1775     if(fV3) fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
1776     if(fV3) fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
1777
1778     // Fill histograms needed for resolution evaluation
1779     if(fV2) fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
1780     if(fV2) fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
1781     if(fV2) fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
1782
1783     if(fV3) fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
1784     if(fV3) fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
1785     if(fV3) fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
1786   }
1787
1788
1789
1790   // clean track array
1791   for(Int_t i=0;i < nusedForK0s;i++){
1792     usedForK0s1[i] = NULL;
1793     usedForK0s2[i] = NULL;
1794   }
1795 }
1796
1797 //_____________________________________________________________________________
1798 Float_t AliAnalysisTaskVnV0::GetVertex(AliAODEvent* aod) const
1799 {
1800
1801   Float_t zvtx = -999;
1802
1803   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
1804   if (!vtxAOD)
1805     return zvtx;
1806   if(vtxAOD->GetNContributors()>0)
1807     zvtx = vtxAOD->GetZ();
1808
1809   return zvtx;
1810 }
1811 //_____________________________________________________________________________
1812 void AliAnalysisTaskVnV0::Terminate(Option_t *)
1813 {
1814   // Terminate loop
1815   Printf("Terminate()");
1816 }
1817 //_____________________________________________________________________________
1818 void AliAnalysisTaskVnV0::OpenInfoCalbration(Int_t run){
1819     TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
1820     TFile *foadb = TFile::Open(oadbfilename.Data());
1821
1822     if(!foadb){
1823         printf("OADB file %s cannot be opened\n",oadbfilename.Data());
1824         return;
1825     }
1826
1827     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1828     if(!cont){
1829         printf("OADB object hMultV0BefCorr is not available in the file\n");
1830         return;
1831     }
1832
1833     if(!(cont->GetObject(run))){
1834         printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
1835         run = 137366;
1836     }
1837     fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
1838
1839     TF1 *fpol0 = new TF1("fpol0","pol0");
1840     fMultV0->Fit(fpol0,"","",0,31);
1841     fV0Cpol = fpol0->GetParameter(0);
1842     fMultV0->Fit(fpol0,"","",32,64);
1843     fV0Apol = fpol0->GetParameter(0);
1844
1845     for(Int_t iside=0;iside<2;iside++){
1846         for(Int_t icoord=0;icoord<2;icoord++){
1847             for(Int_t i=0;i  < 9;i++){
1848                 char namecont[100];
1849                 if(iside==0 && icoord==0)
1850                   snprintf(namecont,100,"hQxc2_%i",i);
1851                 else if(iside==1 && icoord==0)
1852                   snprintf(namecont,100,"hQxa2_%i",i);
1853                 else if(iside==0 && icoord==1)
1854                   snprintf(namecont,100,"hQyc2_%i",i);
1855                 else if(iside==1 && icoord==1)
1856                   snprintf(namecont,100,"hQya2_%i",i);
1857
1858                 cont = (AliOADBContainer*) foadb->Get(namecont);
1859                 if(!cont){
1860                     printf("OADB object %s is not available in the file\n",namecont);
1861                     return;
1862                 }
1863
1864                 if(!(cont->GetObject(run))){
1865                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1866                     run = 137366;
1867                 }
1868                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1869                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1870
1871                 //for v3
1872                 if(iside==0 && icoord==0)
1873                   snprintf(namecont,100,"hQxc3_%i",i);
1874                 else if(iside==1 && icoord==0)
1875                   snprintf(namecont,100,"hQxa3_%i",i);
1876                 else if(iside==0 && icoord==1)
1877                   snprintf(namecont,100,"hQyc3_%i",i);
1878                 else if(iside==1 && icoord==1)
1879                   snprintf(namecont,100,"hQya3_%i",i);
1880
1881                 cont = (AliOADBContainer*) foadb->Get(namecont);
1882                 if(!cont){
1883                     printf("OADB object %s is not available in the file\n",namecont);
1884                     return;
1885                 }
1886
1887                 if(!(cont->GetObject(run))){
1888                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1889                     run = 137366;
1890                 }
1891                 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1892                 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1893
1894             }
1895         }
1896     }
1897 }
1898 //=======================================================================
1899 Int_t AliAnalysisTaskVnV0::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1900 {
1901   Int_t set = 0;
1902   Float_t fV0Cuts[9];
1903   // defines cuts to be used
1904   // fV0Cuts[9] dl dca ctp d0 d0d0 qt minEta maxEta PID
1905   switch(set) {
1906   case(0): // No cuts
1907     fV0Cuts[0] = -1e+6; fV0Cuts[1] = +1e+6; fV0Cuts[2] = -1e+6;
1908     fV0Cuts[3] = -1e+6; fV0Cuts[4] = +1e+6; fV0Cuts[5] = -1e+6;
1909     fV0Cuts[6] = -1e+6; fV0Cuts[7] = +1e+6; fV0Cuts[8] = 0;
1910     break;
1911   case(1): // Tight cuts
1912     fV0Cuts[0] = +0.5; fV0Cuts[1] = +0.5; fV0Cuts[2] = +0.998;
1913     fV0Cuts[3] = +0.1; fV0Cuts[4] = +0.0; fV0Cuts[5] = +0.105;
1914     fV0Cuts[6] = -0.8; fV0Cuts[7] = +0.8; fV0Cuts[8] = 0;
1915     break;
1916   case(2): // Tight cuts + PID
1917     fV0Cuts[0] = +0.5; fV0Cuts[1] = +0.5; fV0Cuts[2] = +0.998;
1918     fV0Cuts[3] = +0.1; fV0Cuts[4] = +0.0; fV0Cuts[5] = +0.105;
1919     fV0Cuts[6] = -0.8; fV0Cuts[7] = +0.8; fV0Cuts[8] = 1;
1920     break;
1921   case(3): // No cuts + PID
1922     fV0Cuts[0] = -1e+6; fV0Cuts[1] = +1e+6; fV0Cuts[2] = -1e+6;
1923     fV0Cuts[3] = -1e+6; fV0Cuts[4] = +1e+6; fV0Cuts[5] = -1e+6;
1924     fV0Cuts[6] = -1e+6; fV0Cuts[7] = +1e+6; fV0Cuts[8] = 1;
1925     break;
1926   }
1927
1928   // daughter cuts
1929   if(! fCutsDaughter){
1930     fCutsDaughter = new AliESDtrackCuts(Form("daughter_cuts_%s","ESD") );
1931     fCutsDaughter->SetPtRange(0.2,10.0);
1932     fCutsDaughter->SetEtaRange(-0.8, 0.8 );
1933     fCutsDaughter->SetMinNClustersTPC(80);
1934     fCutsDaughter->SetMaxChi2PerClusterTPC(4.0);
1935     fCutsDaughter->SetRequireTPCRefit(kTRUE);
1936     fCutsDaughter->SetAcceptKinkDaughters(kFALSE);
1937   }
1938
1939   if (myV0->GetOnFlyStatus() ) return 0;
1940   //the following is needed in order to evualuate track-quality
1941   AliAODTrack *iT, *jT;
1942   AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1943   Double_t pos[3],cov[6];
1944   vV0s->GetXYZ(pos);
1945   vV0s->GetCovarianceMatrix(cov);
1946   const AliESDVertex vESD(pos,cov,100.,100);
1947   // TESTING CHARGE
1948   int iPos, iNeg;
1949   iT=(AliAODTrack*) myV0->GetDaughter(0);
1950   if(iT->Charge()>0) {
1951     iPos = 0; iNeg = 1;
1952   } else {
1953     iPos = 1; iNeg = 0;
1954   }
1955   // END OF TEST
1956
1957   iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1958   AliESDtrack ieT( iT );
1959   ieT.SetTPCClusterMap( iT->GetTPCClusterMap() );
1960   ieT.SetTPCSharedMap( iT->GetTPCSharedMap() );
1961   ieT.SetTPCPointsF( iT->GetTPCNclsF() );
1962   ieT.RelateToVertex(&vESD, tAOD->GetMagneticField(), 100);
1963   if (!fCutsDaughter->IsSelected( &ieT ) ) return 0;
1964
1965   jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1966   AliESDtrack jeT( jT );
1967   jeT.SetTPCClusterMap( jT->GetTPCClusterMap() );
1968   jeT.SetTPCSharedMap( jT->GetTPCSharedMap() );
1969   jeT.SetTPCPointsF( jT->GetTPCNclsF() );
1970   jeT.RelateToVertex(&vESD, tAOD->GetMagneticField(), 100);
1971   if (!fCutsDaughter->IsSelected( &jeT ) ) return 0;
1972
1973   Double_t pvertex[3];
1974   pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1975   pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1976   pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1977   Double_t dDL=myV0->DecayLengthV0( pvertex );
1978   Double_t dDCA=myV0->DcaV0Daughters();
1979   Double_t dCTP=myV0->CosPointingAngle( pvertex );
1980   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1981   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1982   Double_t dD0D0=dD0P*dD0M;
1983   Double_t dQT=myV0->PtArmV0();
1984   Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1985   if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1986 //   Double_t dPT=myV0->Pt();
1987   Double_t dETA=myV0->Eta();
1988   Int_t passes = 1;
1989   if(dDL  <fV0Cuts[0]) passes = 0;
1990   if(dDCA >fV0Cuts[1]) passes = 0;
1991   if(dCTP <fV0Cuts[2]) passes = 0;
1992   if(TMath::Abs(dD0P) <fV0Cuts[3]) passes = 0;
1993   if(TMath::Abs(dD0M) <fV0Cuts[3]) passes = 0;
1994   if(dD0D0>fV0Cuts[4]) passes = 0;
1995   if(dETA <fV0Cuts[6]) passes = 0;
1996   if(dETA >fV0Cuts[7]) passes = 0;
1997   if(specie==0) if(dQT<fV0Cuts[5]) passes = 0;
1998   if(specie==1&&passes==1&&dALPHA<0) passes = 2; // antilambda
1999
2000
2001 //   if(jT->Pt() < 0.5*myV0->Pt() || iT->Pt() < 0.5*myV0->Pt()) passes = 0;
2002
2003
2004   // additional cut
2005 //   if(!(iT->GetStatus() & AliAODTrack::kTPCrefit)) passes = 0;
2006 //   if(!(jT->GetStatus() & AliAODTrack::kTPCrefit)) passes = 0;
2007
2008 //   if(!(iT->GetStatus() & AliAODTrack::kITSrefit)) passes = 0;
2009 //   if(!(jT->GetStatus() & AliAODTrack::kITSrefit)) passes = 0;
2010
2011 //  if(!(iT->GetStatus() & AliAODTrack::kTOFout)) passes = 0;
2012 //  if(!(jT->GetStatus() & AliAODTrack::kTOFout)) passes = 0;
2013
2014   Bool_t trkFlag = iT->TestFilterBit(1); // TPC only tracks (4,global track)
2015   Bool_t trkFlag2 = jT->TestFilterBit(1); // TPC only tracks (4,global track)
2016
2017   if(!trkFlag) passes = 0;
2018   if(!trkFlag2) passes = 0;
2019
2020   if(passes&&fV0Cuts[8]) {
2021
2022     Double_t dedxExp[8];
2023     fPID->ComputeProb(iT,tAOD); // compute Bayesian probabilities
2024     Float_t nsigmaTPC[8];
2025
2026     Int_t tofMatch=0;
2027     Int_t tofMatch2=0;
2028
2029     if(iT->GetDetPid()){ // check the PID object is available
2030       for(Int_t iS=0;iS < 8;iS++){
2031         dedxExp[iS] = fPID->GetExpDeDx(iT,iS);
2032         nsigmaTPC[iS] = (fPID->GetDeDx() - dedxExp[iS])/(dedxExp[iS]*0.07);
2033       }
2034     }
2035     else{
2036       for(Int_t iS=0;iS < 8;iS++)
2037         nsigmaTPC[iS] = 10;
2038     }
2039
2040
2041     if(fPID->GetCurrentMask(1)) // if TOF is present
2042       tofMatch = 1;
2043
2044 //     Float_t tofMismProbMC = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2045
2046     Float_t *probRead = fPID->GetProb();
2047     Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2048
2049     fPID->ComputeProb(jT,tAOD); // compute Bayesian probabilities
2050     Float_t nsigmaTPC2[8];
2051     if(jT->GetDetPid()){ // check the PID object is available
2052       for(Int_t iS=0;iS < 8;iS++){
2053         dedxExp[iS] = fPID->GetExpDeDx(jT,iS);
2054         nsigmaTPC2[iS] = (fPID->GetDeDx() - dedxExp[iS])/(dedxExp[iS]*0.07);
2055       }
2056     }
2057     else{
2058       for(Int_t iS=0;iS < 8;iS++)
2059         nsigmaTPC2[iS] = 10;
2060     }
2061
2062     if(fPID->GetCurrentMask(1)) // if TOF is present
2063       tofMatch2 = 1;
2064
2065 //     Float_t tofMismProbMC2 = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2066
2067     probRead = fPID->GetProb();
2068     Float_t prob2[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2069
2070     if(jT->GetTPCNcls() < fNcluster) passes = 0;
2071     else if(iT->GetTPCNcls() < fNcluster) passes = 0;
2072
2073 //     if(! (tofMatch && tofMatch2)) passes = 0;
2074
2075     /*
2076     Float_t dMASS = myV0->MassK0Short();
2077     Float_t nsigmaMass = TMath::Abs(dMASS-0.497)/0.005;
2078     if(specie == 0 && TMath::Abs(nsigmaMass) < 1 && myV0->Pt() > 1) printf("candidate i=(pt=%f-phi=%f-tof=%i) j=(pt=%f-phi=%f-tof=%i) \n",iT->Pt(),iT->Phi(),tofMatch,jT->Pt(),jT->Phi(),tofMatch2);
2079     */
2080
2081     switch(specie) {
2082     case 0: // K0 PID
2083       if(0){
2084         if( ((jT->GetTPCmomentum()<15) &&
2085              (TMath::Abs(nsigmaTPC2[2])>3.)) || prob2[2] < 0.9)
2086           passes = 0;
2087         if( ((iT->GetTPCmomentum()<15) &&
2088              (TMath::Abs(nsigmaTPC[2])>3.))|| prob[2] < 0.9 )
2089           passes = 0;
2090       }
2091       break;
2092     case 1: // Lambda PID  i==pos j ==neg
2093       if(passes==1) {
2094         if( (iT->GetTPCmomentum()<15) &&
2095             (TMath::Abs(nsigmaTPC[4])>3.) )
2096           passes = 0;
2097         if( (jT->GetTPCmomentum()<15) &&
2098             (TMath::Abs(nsigmaTPC2[2])>3.) )
2099           passes = 0;
2100       }
2101       if(passes==2) {
2102         if( (iT->GetTPCmomentum()<15) &&
2103             (TMath::Abs(nsigmaTPC[2])>3.) )
2104           passes = 0;
2105         if( (jT->GetTPCmomentum()<15) &&
2106             (TMath::Abs(nsigmaTPC2[4])>3.) )
2107           passes = 0;
2108       }
2109       break;
2110     }
2111   }
2112   return passes;
2113 }
2114 //=======================================================================
2115 void AliAnalysisTaskVnV0::SelectK0s(){
2116   fNK0s=0;
2117   fNpiPos=0;
2118   fNpiNeg=0;
2119
2120   if(fModulationDEDx) fPID->SetPsiCorrectionDeDx(evPlAng2,1.0); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed)
2121
2122   // fill pion stacks
2123   Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
2124   for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
2125     AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
2126
2127     if (!aodTrack){
2128       aodTrack->Delete();
2129       continue;
2130     }
2131
2132     Bool_t trkFlag = aodTrack->TestFilterBit(1); // TPC only tracks
2133 //    trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
2134
2135     if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
2136       continue;
2137     }
2138
2139     Double_t b[2] = {-99., -99.};
2140     Double_t bCov[3] = {-99., -99., -99.};
2141     if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
2142       continue;
2143
2144     if(TMath::Abs(b[0]) < 0.5/aodTrack->Pt()) continue;
2145
2146     fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
2147     Float_t *probRead = fPID->GetProb();
2148     Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2149  //    Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2150
2151     Int_t charge = aodTrack->Charge();
2152     if(prob[2] > 0.9){
2153       if(charge > 0){
2154         fIPiPos[fNpiPos] = iT;
2155         fNpiPos++;
2156       }
2157       else{
2158         fIPiNeg[fNpiNeg] = iT;
2159         fNpiNeg++;
2160       }
2161     }
2162   }
2163
2164   for(Int_t i=0;i < fNpiPos;i++){
2165     AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
2166     AliESDtrack pipE(pip);
2167
2168     for(Int_t j=0;j < fNpiNeg;j++){
2169       AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
2170       AliESDtrack pinE(pin);
2171
2172       Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
2173
2174       Double_t pPos[3];
2175       Double_t pNeg[3];
2176       pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
2177       pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
2178
2179       Float_t length = (xp+xn)*0.5;
2180
2181       Float_t pxs = pPos[0] + pNeg[0];
2182       Float_t pys = pPos[1] + pNeg[1];
2183       Float_t pzs = pPos[2] + pNeg[2];
2184       Float_t es = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.13957*0.13957) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.13957*0.13957);
2185
2186       Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
2187       Float_t phi = TMath::ATan2(pys,pxs);
2188       Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
2189
2190       //      if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
2191
2192       if(mindist < 0.2&& length > 1 && length < 25){
2193         fHK0sMass->Fill(pt,mass);
2194
2195         Float_t esL = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.938*0.938) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.13957*0.13957);
2196         Float_t esAL = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.13957*0.13957) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.938*0.938);
2197
2198         Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
2199         Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
2200
2201         fHK0vsLambda->Fill(mass,TMath::Min(massaL,massaAL));
2202
2203         if(TMath::Abs(mass-0.497)/0.005 < 1 && massaL > 1.15 && massaAL > 1.15){
2204           fPhiK0s[fNK0s] = phi;
2205           fPtK0s[fNK0s] = pt;
2206           fNK0s++;
2207         }
2208       }
2209     }
2210   }
2211 }