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