From afa8df58bfb681136374c86f82e07721a5840397 Mon Sep 17 00:00:00 2001 From: iseliouj Date: Mon, 27 Feb 2012 14:19:33 +0000 Subject: [PATCH] From Francesco: Added AliAnalysisTaskVnV0 and macros --- PWG/CMakelibPWGflowTasks.pkg | 1 + PWG/FLOW/Tasks/AliAnalysisTaskVnV0.cxx | 1062 ++++++++++++++++++++++++ PWG/FLOW/Tasks/AliAnalysisTaskVnV0.h | 102 +++ PWG/PWGflowTasksLinkDef.h | 1 + PWGCF/FLOW/macros/AddTaskVZERO.C | 45 + PWGCF/FLOW/macros/extractFlowVZERO.C | 279 +++++++ 6 files changed, 1490 insertions(+) create mode 100644 PWG/FLOW/Tasks/AliAnalysisTaskVnV0.cxx create mode 100644 PWG/FLOW/Tasks/AliAnalysisTaskVnV0.h create mode 100644 PWGCF/FLOW/macros/AddTaskVZERO.C create mode 100644 PWGCF/FLOW/macros/extractFlowVZERO.C diff --git a/PWG/CMakelibPWGflowTasks.pkg b/PWG/CMakelibPWGflowTasks.pkg index f35d70063c5..176f139aa7f 100644 --- a/PWG/CMakelibPWGflowTasks.pkg +++ b/PWG/CMakelibPWGflowTasks.pkg @@ -51,6 +51,7 @@ set ( SRCS FLOW/Tasks/AliFlowBayesianPID.cxx FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx FLOW/Tasks/AliAnalysisTaskFilterFE.cxx + FLOW/Tasks/AliAnalysisTaskVnV0.cxx ) string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" ) diff --git a/PWG/FLOW/Tasks/AliAnalysisTaskVnV0.cxx b/PWG/FLOW/Tasks/AliAnalysisTaskVnV0.cxx new file mode 100644 index 00000000000..10c018f8f2e --- /dev/null +++ b/PWG/FLOW/Tasks/AliAnalysisTaskVnV0.cxx @@ -0,0 +1,1062 @@ +#include "AliAnalysisTaskVnV0.h" + +// ROOT includes +#include + +// AliRoot includes +#include "AliInputEventHandler.h" +#include "AliAODEvent.h" +#include "AliAODVertex.h" +#include "AliAODTrack.h" +#include "AliCentrality.h" +#include "AliVHeader.h" +#include "AliAODVZERO.h" +#include "TFile.h" +#include "AliOADBContainer.h" +#include "TH2F.h" +#include "TF1.h" + +// STL includes +//#include +//using namespace std; + +ClassImp(AliAnalysisTaskVnV0) + +/* +fMeanQ[nCentrBin][2][2]; // and recentering + Float_t fWidthQ[nCentrBin][2][2]; // ... + Float_t fMeanQv3[nCentrBin][2][2]; // also for v3 + Float_t fWidthQv3[nCentrBin][2][2]; // ... +*/ + +//_____________________________________________________________________________ +AliAnalysisTaskVnV0::AliAnalysisTaskVnV0(): + AliAnalysisTaskSE(), + fAOD(0), + fVtxCut(10), // cut on |vertex| < fVtxCut + fEtaCut(0.8), // cut on |eta| < fEtaCut + fMinPt(0.15), // cut on pt > fMinPt + fRun(-1), + fList(new TList()), + fList2(new TList()), + fMultV0(0), + fV0Cpol(100), + fV0Apol(100), + fContAllChargesV0A(0), + fContAllChargesV0C(0), + fContAllChargesV0Av3(0), + fContAllChargesV0Cv3(0), + fHResTPCv0A2(0), + fHResTPCv0C2(0), + fHResv0Cv0A2(0), + fHResTPCv0A3(0), + fHResTPCv0C3(0), + fHResv0Cv0A3(0), + fPhiRPv0A(0), + fPhiRPv0C(0), + fPhiRPv0Av3(0), + fPhiRPv0Cv3(0), + fPhiTracks(0), + fQA(0), + fQA2(0), + fQAv3(0), + fQA2v3(0), + fPID(new AliFlowBayesianPID()), + fTree(0), + fCentrality(0), + evPlAngV0ACor2(0), + evPlAngV0CCor2(0), + evPlAng2(0), + evPlAngV0ACor3(0), + evPlAngV0CCor3(0), + evPlAng3(0), + fV2(kTRUE), + fV3(kTRUE) +{ + DefineOutput(1, TList::Class()); + DefineOutput(2, TList::Class()); + + // Default constructor (should not be used) + fList->SetName("resultsV2"); + fList2->SetName("resultsV3"); + + fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h + +} + +//______________________________________________________________________________ +AliAnalysisTaskVnV0::AliAnalysisTaskVnV0(const char *name): + AliAnalysisTaskSE(name), + fAOD(0), + fVtxCut(10), // cut on |vertex| < fVtxCut + fEtaCut(0.8), // cut on |eta| < fEtaCut + fMinPt(0.15), // cut on pt > fMinPt + fRun(-1), + fList(new TList()), + fList2(new TList()), + fMultV0(0), + fV0Cpol(100), + fV0Apol(100), + fContAllChargesV0A(0), + fContAllChargesV0C(0), + fContAllChargesV0Av3(0), + fContAllChargesV0Cv3(0), + fHResTPCv0A2(0), + fHResTPCv0C2(0), + fHResv0Cv0A2(0), + fHResTPCv0A3(0), + fHResTPCv0C3(0), + fHResv0Cv0A3(0), + fPhiRPv0A(0), + fPhiRPv0C(0), + fPhiRPv0Av3(0), + fPhiRPv0Cv3(0), + fPhiTracks(0), + fQA(0), + fQA2(0), + fQAv3(0), + fQA2v3(0), + fPID(new AliFlowBayesianPID()), + fTree(0), + fCentrality(0), + evPlAngV0ACor2(0), + evPlAngV0CCor2(0), + evPlAng2(0), + evPlAngV0ACor3(0), + evPlAngV0CCor3(0), + evPlAng3(0), + fV2(kTRUE), + fV3(kTRUE) +{ + + DefineOutput(1, TList::Class()); + DefineOutput(2, TList::Class()); + + // Output slot #1 writes into a TTree + fList->SetName("resultsV2"); + fList2->SetName("resultsV3"); + + fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h +} + +//_____________________________________________________________________________ +AliAnalysisTaskVnV0::~AliAnalysisTaskVnV0() +{ + +} + +//______________________________________________________________________________ +void AliAnalysisTaskVnV0::UserCreateOutputObjects() +{ + + // Tree for EP debug (comment the adding to v2 list id not needed) + fTree = new TTree("tree","tree"); + fTree->Branch("evPlAngV0ACor2",&evPlAngV0ACor2,"evPlAngV0ACor2/F"); + fTree->Branch("evPlAngV0CCor2",&evPlAngV0CCor2,"evPlAngV0CCor2/F"); + fTree->Branch("evPlAng2",&evPlAng2,"evPlAng2/F"); + fTree->Branch("fCentrality",&fCentrality,"fCentrality/F"); + fTree->Branch("evPlAngV0ACor3",&evPlAngV0ACor3,"evPlAngV0ACor3/F"); + fTree->Branch("evPlAngV0CCor3",&evPlAngV0CCor3,"evPlAngV0CCor3/F"); + fTree->Branch("evPlAng3",&evPlAng3,"evPlAng3/F"); + + + // Container analyses (different steps mean different species) + const Int_t nPtBinsTOF = 45; + 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}; + const Int_t nChargeBinsTOF = 3; + Double_t binChargeTOF[nChargeBinsTOF+1] = {-1.5,-0.5,0.5,1.5}; + const Int_t nv_2BinsTOF = 50; + Double_t binV_2TOF[nv_2BinsTOF+1]; + for(Int_t i=0;iSetBinLimits(0,binCentrTOF); + fContAllChargesV0A->SetBinLimits(1,binV_2TOF); + fContAllChargesV0A->SetBinLimits(2,binChargeTOF); + fContAllChargesV0A->SetBinLimits(3,binsPtTOF); + fContAllChargesV0A->SetBinLimits(4,binProbTOF); + fContAllChargesV0A->SetBinLimits(5,binPsiTOF); + fContAllChargesV0A->SetBinLimits(6,binMaskPID); + fContAllChargesV0C->SetBinLimits(0,binCentrTOF); + fContAllChargesV0C->SetBinLimits(1,binV_2TOF); + fContAllChargesV0C->SetBinLimits(2,binChargeTOF); + fContAllChargesV0C->SetBinLimits(3,binsPtTOF); + fContAllChargesV0C->SetBinLimits(4,binProbTOF); + fContAllChargesV0C->SetBinLimits(5,binPsiTOF); + fContAllChargesV0C->SetBinLimits(6,binMaskPID); + + fContAllChargesV0A->SetVarTitle(0,"centrality"); + fContAllChargesV0A->SetVarTitle(1,"v_{2}"); + fContAllChargesV0A->SetVarTitle(2,"charge"); + fContAllChargesV0A->SetVarTitle(3,"p_{t} (GeV/c)"); + fContAllChargesV0A->SetVarTitle(4,"Bayesian Probability"); + fContAllChargesV0A->SetVarTitle(5,"#Psi V0A"); + fContAllChargesV0A->SetVarTitle(6,"TOF PID"); + fContAllChargesV0C->SetVarTitle(0,"centrality"); + fContAllChargesV0C->SetVarTitle(1,"v_{2}"); + fContAllChargesV0C->SetVarTitle(2,"charge"); + fContAllChargesV0C->SetVarTitle(3,"p_{t} (GeV/c)"); + fContAllChargesV0C->SetVarTitle(4,"Bayesian Probability"); + fContAllChargesV0C->SetVarTitle(5,"#Psi V0C"); + fContAllChargesV0C->SetVarTitle(6,"TOF PID"); + + fList->Add(fContAllChargesV0A); + fList->Add(fContAllChargesV0C); + + // v3 container + fContAllChargesV0Av3 = new AliCFContainer("fContAllChargesV0Av3","centr:v3:charge:pt:prob:Psi:PIDmask",1+7,7,binsTOF); + fContAllChargesV0Cv3 = new AliCFContainer("fContAllChargesV0Cv3","centr:v3:charge:pt:prob:Psi:PIDmask",1+7,7,binsTOF); + fContAllChargesV0Av3->SetBinLimits(0,binCentrTOF); + fContAllChargesV0Av3->SetBinLimits(1,binV_2TOF); + fContAllChargesV0Av3->SetBinLimits(2,binChargeTOF); + fContAllChargesV0Av3->SetBinLimits(3,binsPtTOF); + fContAllChargesV0Av3->SetBinLimits(4,binProbTOF); + fContAllChargesV0Av3->SetBinLimits(5,binPsiTOFv3); + fContAllChargesV0Av3->SetBinLimits(6,binMaskPID); + fContAllChargesV0Cv3->SetBinLimits(0,binCentrTOF); + fContAllChargesV0Cv3->SetBinLimits(1,binV_2TOF); + fContAllChargesV0Cv3->SetBinLimits(2,binChargeTOF); + fContAllChargesV0Cv3->SetBinLimits(3,binsPtTOF); + fContAllChargesV0Cv3->SetBinLimits(4,binProbTOF); + fContAllChargesV0Cv3->SetBinLimits(5,binPsiTOFv3); + fContAllChargesV0Cv3->SetBinLimits(6,binMaskPID); + + fContAllChargesV0Av3->SetVarTitle(0,"centrality"); + fContAllChargesV0Av3->SetVarTitle(1,"v_{2}"); + fContAllChargesV0Av3->SetVarTitle(2,"charge"); + fContAllChargesV0Av3->SetVarTitle(3,"p_{t} (GeV/c)"); + fContAllChargesV0Av3->SetVarTitle(4,"Bayesian Probability"); + fContAllChargesV0Av3->SetVarTitle(5,"#Psi V0A"); + fContAllChargesV0Av3->SetVarTitle(6,"TOF PID"); + fContAllChargesV0Cv3->SetVarTitle(0,"centrality"); + fContAllChargesV0Cv3->SetVarTitle(1,"v_{3}"); + fContAllChargesV0Cv3->SetVarTitle(2,"charge"); + fContAllChargesV0Cv3->SetVarTitle(3,"p_{t} (GeV/c)"); + fContAllChargesV0Cv3->SetVarTitle(4,"Bayesian Probability"); + fContAllChargesV0Cv3->SetVarTitle(5,"#Psi V0C"); + fContAllChargesV0Cv3->SetVarTitle(6,"TOF PID"); + + fList2->Add(fContAllChargesV0Av3); + fList2->Add(fContAllChargesV0Cv3); + + // TProfile for resolutions 3 subevents (V0A, V0C, TPC) + // v2 + fHResTPCv0A2 = new TProfile("hResTPCv0A2","",9,0,9); + fHResTPCv0C2 = new TProfile("hResTPCv0C2","",9,0,9); + fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",9,0,9); + + fList->Add(fHResTPCv0A2); + fList->Add(fHResTPCv0C2); + fList->Add(fHResv0Cv0A2); + + // v3 + fHResTPCv0A3 = new TProfile("hResTPCv0A3","",9,0,9); + fHResTPCv0C3 = new TProfile("hResTPCv0C3","",9,0,9); + fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",9,0,9); + + fList2->Add(fHResTPCv0A3); + fList2->Add(fHResTPCv0C3); + fList2->Add(fHResv0Cv0A3); + + // V0A and V0C event plane distributions + //v2 + fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi(),TMath::Pi()); + fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi(),TMath::Pi()); + + //v3 + fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi()/3*2,TMath::Pi()/3*2); + fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi()/3*2,TMath::Pi()/3*2); + + // Track phi distribution container (only once) put in v2 output (only if needed) + const Int_t nPhiBin = 20; + Double_t binPhi[nPhiBin+1]; + for(Int_t i=0;iSetBinLimits(0,binCentrTOF); + fPhiTracks->SetBinLimits(1,binChargeTOF); + fPhiTracks->SetBinLimits(2,binsPtTOF); + fPhiTracks->SetBinLimits(3,binProbTOF); + fPhiTracks->SetBinLimits(4,binPhi); + fPhiTracks->SetVarTitle(0,"centrality"); + fPhiTracks->SetVarTitle(1,"charge"); + fPhiTracks->SetVarTitle(2,"p_{t} (GeV/c)"); + fPhiTracks->SetVarTitle(3,"Bayesian Probability"); + fPhiTracks->SetVarTitle(4,"#phi"); + + // QA container + // v2 + const Int_t nDETsignal = 50; + Double_t binDETsignal[nDETsignal+1]; + for(Int_t i=0;iSetBinLimits(0,binCentrTOF); + fQA->SetBinLimits(1,binsPtTOF); + fQA->SetBinLimits(2,binProbTOF); + fQA->SetBinLimits(3,binDETsignal); + fQA->SetBinLimits(4,binDETsignal); + fQA->SetBinLimits(5,binDeltaPhi); + fQA->SetBinLimits(6,binMaskPID); + fQA->SetVarTitle(0,"centrality"); + fQA->SetVarTitle(1,"p_{t} (GeV/c)"); + fQA->SetVarTitle(2,"Bayesian Probability"); + fQA->SetVarTitle(3,"N_{#sigma}^{TPC}"); + fQA->SetVarTitle(4,"N_{#sigma}^{TOF}"); + fQA->SetVarTitle(5,"#Delta#phi (V0A)"); + fQA->SetVarTitle(6,"TOF PID"); + + fQA2 = new AliCFContainer("fQA2v2","centr:pt:prob:TPCsig:TOFsig:DeltaPhi:maskPID",7,8,binsQA); + fQA2->SetBinLimits(0,binCentrTOF); + fQA2->SetBinLimits(1,binsPtTOF); + fQA2->SetBinLimits(2,binProbTOF); + fQA2->SetBinLimits(3,binDETsignal); + fQA2->SetBinLimits(4,binDETsignal); + fQA2->SetBinLimits(5,binDeltaPhi); + fQA2->SetBinLimits(6,binMaskPID); + fQA2->SetVarTitle(0,"centrality"); + fQA2->SetVarTitle(1,"p_{t} (GeV/c)"); + fQA2->SetVarTitle(2,"Bayesian Probability"); + fQA2->SetVarTitle(3,"N_{#sigma}^{TPC}"); + fQA2->SetVarTitle(4,"N_{#sigma}^{TOF}"); + fQA2->SetVarTitle(5,"#Delta#phi (V0C)"); + fQA2->SetVarTitle(6,"TOF PID"); + + // v3 + const Int_t nDeltaPhiV3 = 7; + Double_t binDeltaPhiV3[nDeltaPhiV3+1]; + for(Int_t i=0;iSetBinLimits(0,binCentrTOF); + fQAv3->SetBinLimits(1,binsPtTOF); + fQAv3->SetBinLimits(2,binProbTOF); + fQAv3->SetBinLimits(3,binDETsignal); + fQAv3->SetBinLimits(4,binDETsignal); + fQAv3->SetBinLimits(5,binDeltaPhiV3); + fQAv3->SetBinLimits(6,binMaskPID); + fQAv3->SetVarTitle(0,"centrality"); + fQAv3->SetVarTitle(1,"p_{t} (GeV/c)"); + fQAv3->SetVarTitle(2,"Bayesian Probability"); + fQAv3->SetVarTitle(3,"N_{#sigma}^{TPC}"); + fQAv3->SetVarTitle(4,"N_{#sigma}^{TOF}"); + fQAv3->SetVarTitle(5,"#Delta#phi (V0A)"); + fQAv3->SetVarTitle(6,"TOF PID"); + + fQA2v3 = new AliCFContainer("fQA2v3","centr:pt:prob:TPCsig:TOFsig:DeltaPhi:maskPID",7,8,binsQA); + fQA2v3->SetBinLimits(0,binCentrTOF); + fQA2v3->SetBinLimits(1,binsPtTOF); + fQA2v3->SetBinLimits(2,binProbTOF); + fQA2v3->SetBinLimits(3,binDETsignal); + fQA2v3->SetBinLimits(4,binDETsignal); + fQA2v3->SetBinLimits(5,binDeltaPhiV3); + fQA2v3->SetBinLimits(6,binMaskPID); + fQA2v3->SetVarTitle(0,"centrality"); + fQA2v3->SetVarTitle(1,"p_{t} (GeV/c)"); + fQA2v3->SetVarTitle(2,"Bayesian Probability"); + fQA2v3->SetVarTitle(3,"N_{#sigma}^{TPC}"); + fQA2v3->SetVarTitle(4,"N_{#sigma}^{TOF}"); + fQA2v3->SetVarTitle(5,"#Delta#phi (V0C)"); + fQA2v3->SetVarTitle(6,"TOF PID"); + + + fList->Add(fPhiRPv0A); + fList->Add(fPhiRPv0C); + // fList->Add(fPhiTracks); // comment if not needed + fList->Add(fQA); + fList->Add(fQA2); + + fList2->Add(fPhiRPv0Av3); + fList2->Add(fPhiRPv0Cv3); + fList2->Add(fQAv3); + fList2->Add(fQA2v3); + + // fList->Add(fTree); // comment if not needed + + printf("Output creation ok!!\n\n\n\n"); + + // Post output data. + if(fV2) PostData(1, fList); + if(fV3) PostData(2, fList2); +} + +//______________________________________________________________________________ +void AliAnalysisTaskVnV0::UserExec(Option_t *) +{ + // Main loop + // Called for each event + + fAOD = dynamic_cast(InputEvent()); + if(!fAOD){ + Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__); + this->Dump(); + return; + } + + Int_t run = fAOD->GetRunNumber(); + + if(run != fRun){ + // Load the calibrations run dependent + OpenInfoCalbration(run); + fRun=run; + } + + Float_t zvtx = GetVertex(fAOD); + + if (TMath::Abs(zvtx) < fVtxCut) { + //Centrality + Float_t v0Centr = -10.; + Float_t trkCentr = -10.; + AliCentrality *centrality = fAOD->GetCentrality(); + if (centrality){ + v0Centr = centrality->GetCentralityPercentile("V0M"); + trkCentr = centrality->GetCentralityPercentile("TRK"); + } + + if(TMath::Abs(v0Centr - trkCentr) < 5.0){ // consistency cut on centrality selection + fPID->SetDetResponse(fAOD, v0Centr); // Set the PID object for each event!!!! + Analyze(fAOD,v0Centr); // Do analysis!!!! + + fCentrality = v0Centr; + if(fV2) fTree->Fill(); + } + } + +} + +//________________________________________________________________________ +void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr) +{ + 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}; + + // Event plane resolution for v2 + 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 + 0.446480,0.612705,0.712222,0.736200,0.697907,0.610114,0.481009,0.327402,0.182277};// V0C vs. centrality + + Int_t iC = -1; + if (v0Centr >0 && v0Centr < 80){ // analysis only for 0-80% centrality classes + // centrality bins + if(v0Centr < 5) iC = 0; + else if(v0Centr < 10) iC = 1; + else if(v0Centr < 20) iC = 2; + else if(v0Centr < 30) iC = 3; + else if(v0Centr < 40) iC = 4; + else if(v0Centr < 50) iC = 5; + else if(v0Centr < 60) iC = 6; + else if(v0Centr < 70) iC = 7; + else iC = 8; + + //reset Q vector info + Double_t Qxa2 = 0, Qya2 = 0; + Double_t Qxc2 = 0, Qyc2 = 0; + Double_t Qxa3 = 0, Qya3 = 0; + Double_t Qxc3 = 0, Qyc3 = 0; + + //V0 info + AliAODVZERO* aodV0 = aodEvent->GetVZEROData(); + + for (Int_t iv0 = 0; iv0 < 64; iv0++) { + Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8); + Float_t multv0 = aodV0->GetMultiplicity(iv0); + if (iv0 < 32){ // V0C + Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1); + Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1); + Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1); + Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1); + } else { // V0A + Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1); + Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1); + Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1); + Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1); + } + } + + //grab for each centrality the proper histo with the Qx and Qy to do the recentering + Double_t Qxamean2 = fMeanQ[iC][1][0]; + Double_t Qxarms2 = fWidthQ[iC][1][0]; + Double_t Qyamean2 = fMeanQ[iC][1][1]; + Double_t Qyarms2 = fWidthQ[iC][1][1]; + Double_t Qxamean3 = fMeanQv3[iC][1][0]; + Double_t Qxarms3 = fWidthQv3[iC][1][0]; + Double_t Qyamean3 = fMeanQv3[iC][1][1]; + Double_t Qyarms3 = fWidthQv3[iC][1][1]; + + Double_t Qxcmean2 = fMeanQ[iC][0][0]; + Double_t Qxcrms2 = fWidthQ[iC][0][0]; + Double_t Qycmean2 = fMeanQ[iC][0][1]; + Double_t Qycrms2 = fWidthQ[iC][0][1]; + Double_t Qxcmean3 = fMeanQv3[iC][0][0]; + Double_t Qxcrms3 = fWidthQv3[iC][0][0]; + Double_t Qycmean3 = fMeanQv3[iC][0][1]; + Double_t Qycrms3 = fWidthQv3[iC][0][1]; + + Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2; + Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2; + Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2; + Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2; + Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3; + Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3; + Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3; + Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3; + + evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.; + evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.; + evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.; + evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.; + + Int_t nAODTracks = aodEvent->GetNumberOfTracks(); + //loop track and get pid + for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks + AliAODTrack* aodTrack = aodEvent->GetTrack(iT); + + if (!aodTrack){ + aodTrack->Delete(); + continue; + } + + Bool_t trkFlag = aodTrack->TestFilterBit(1); // TPC only tracks + + if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < 70) || !trkFlag){ + continue; + } + + Double_t b[2] = {-99., -99.}; + Double_t bCov[3] = {-99., -99., -99.}; + if (!aodTrack->PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) + continue; + + if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4)) + continue; + + // re-map the container in an array to do the analysis for V0A and V0C within a loop + Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2}; + AliCFContainer *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C}; + AliCFContainer *QA[2] = {fQA,fQA2}; + + Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3}; + AliCFContainer *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3}; + AliCFContainer *QAv3[2] = {fQAv3,fQA2v3}; + + for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side + + fPID->SetPsiCorrectionDeDx(evPlAngV0[iV0],evPlRes[iV0*8+iC]); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed) + + Double_t v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0])); + Double_t v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0])); + + fPID->ComputeProb(aodTrack); // compute Bayesian probabilities + Float_t *probRead = fPID->GetProb(); + Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]}; + Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis + Double_t x[7] = {iC,v2V0,aodTrack->Charge(),aodTrack->Pt(),1,evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5}; // to fill analysis v2 container + Double_t x3[7] = {iC,v3V0,aodTrack->Charge(),aodTrack->Pt(),1,evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5}; // to fill analysis v3 container + + Double_t phi[5] = {iC,aodTrack->Charge(),aodTrack->Pt(),1,aodTrack->Phi()}; // to fill track container + + // Fill no PID + if(iV0 && fV2) fPhiTracks->Fill(phi,0); + if(fV2) contV0[iV0]->Fill(x,0); + if(fV3) contV0v3[iV0]->Fill(x3,0); + + + Float_t dedx = 0; + Double_t dedxExp[8]; + Float_t tof = -1; + Double_t inttimes[8] = {0.,0.,0.,0.,0.,0.,0.,0.}; + Double_t expTOFsigma[8] = {0.,0.,0.,0.,0.,0.,0.,0.}; + if(aodTrack->GetDetPid()){ // check the PID object is available + for(Int_t iS=0;iS < 8;iS++) + dedxExp[iS] = fPID->GetExpDeDx(aodTrack,iS); + + dedx = aodTrack->GetTPCsignal(); + + if(fPID->GetCurrentMask(1)){ // if TOF is present + Float_t ptrack = aodTrack->P(); + tof = aodTrack->GetTOFsignal() - fPID->GetESDpid()->GetTOFResponse().GetStartTime(ptrack); + aodTrack->GetIntegratedTimes(inttimes); + + for(Int_t iS=5;iS < 8;iS++) // extra info for light nuclei + inttimes[iS] = inttimes[0] / ptrack * mass[iS] * TMath::Sqrt(1+ptrack*ptrack/mass[iS]/mass[iS]); + + for(Int_t iS=0;iS<8;iS++) expTOFsigma[iS] = fPID->GetESDpid()->GetTOFResponse().GetExpectedSigma(ptrack, inttimes[iS], mass[iS]); + } + } + + Float_t deltaPhiV0 = aodTrack->Phi() - evPlAngV0[iV0]; + if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi(); + else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi(); + if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi(); + else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi(); + + Float_t deltaPhiV0v3 = aodTrack->Phi() - evPlAngV0v3[iV0]; + if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi(); + else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi(); + if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi(); + else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi(); + + // variable to fill QA container + Double_t xQA[7] = {iC,aodTrack->Pt(), 0.0, 4.99, 4.99,deltaPhiV0,x[6]}; // v2 + Double_t xQA3[7] = {iC,aodTrack->Pt(), 0.0, 4.99, 4.99,deltaPhiV0v3,x[6]}; // v3 + + //pid selection + if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID and PID object strictly required (very important!!!!) + else if(prob[2] > 0.6){ // pi + phi[3] = prob[2]; // set probability in the container variables + x[4] = prob[2]; + xQA[2] = prob[2]; + x3[4] = x[4]; + xQA3[2] = xQA[2]; + if(dedx > 10.){ // set TPC signal in the QA container variables + xQA[3] = (dedx - dedxExp[2])/(dedxExp[2]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(x[6] > 0.5){ // set TOF signal in the QA container variables if present + xQA[4] = (tof - inttimes[2])/expTOFsigma[2]; // TOF + xQA3[4] = xQA[4]; // TOF + } + if(TMath::Abs(xQA[3]) < 5){ // TPC 5 sigma extra cut to accept the track + if(iV0 && fV2) fPhiTracks->Fill(phi,1); + if(fV2) contV0[iV0]->Fill(x,1); + if(fV2) QA[iV0]->Fill(xQA,0); + if(fV3) contV0v3[iV0]->Fill(x3,1); + if(fV3) QAv3[iV0]->Fill(xQA3,0); + } + } + else if(prob[3] > 0.6){ // K + phi[3] = prob[3]; + x[4] = prob[3]; + xQA[2] = prob[3]; + x3[4] = x[4]; + xQA3[2] = xQA[2]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[3])/(dedxExp[3]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(x[6] > 0.5){ + xQA[4] = (tof - inttimes[3])/expTOFsigma[3]; // TOF + xQA3[4] = xQA[4]; // TOF + } + if(TMath::Abs(xQA[3]) < 5){ + if(iV0 && fV2) fPhiTracks->Fill(phi,2); + if(fV2) contV0[iV0]->Fill(x,2); + if(fV2) QA[iV0]->Fill(xQA,1); + if(fV3) contV0v3[iV0]->Fill(x3,2); + if(fV3) QAv3[iV0]->Fill(xQA3,1); + } + } + else if(prob[4] > 0.6){ // p + phi[3] = prob[4]; + x[4] = prob[4]; + xQA[2] = prob[4]; + x3[4] = x[4]; + xQA3[2] = xQA[2]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[4])/(dedxExp[4]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(x[6] > 0.5){ + xQA[4] = (tof - inttimes[4])/expTOFsigma[4]; // TOF + xQA3[4] = xQA[4]; // TOF + } + if(TMath::Abs(xQA[3]) < 5){ + if(iV0 && fV2) fPhiTracks->Fill(phi,3); + if(fV2) contV0[iV0]->Fill(x,3); + if(fV2) QA[iV0]->Fill(xQA,2); + if(fV3) contV0v3[iV0]->Fill(x3,3); + if(fV3) QAv3[iV0]->Fill(xQA3,2); + } + } + else if(prob[0] > 0.6){ // e + phi[3] = prob[0]; + x[4] = prob[0]; + xQA[2] = prob[0]; + x3[4] = x[4]; + xQA3[2] = xQA[2]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[0])/(dedxExp[0]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(x[6] > 0.5){ + xQA[4] = (tof - inttimes[0])/expTOFsigma[0]; // TOF + xQA3[4] = xQA[4]; // TOF + } + if(TMath::Abs(xQA[3]) < 5){ + if(iV0 && fV2) fPhiTracks->Fill(phi,4); + if(fV2) contV0[iV0]->Fill(x,4); + if(fV2) QA[iV0]->Fill(xQA,3); + if(fV3) contV0v3[iV0]->Fill(x3,4); + if(fV3) QAv3[iV0]->Fill(xQA3,3); + } + } + else if(prob[5] > 0.6){ // d + phi[3] = prob[5]; + x[4] = prob[5]; + xQA[2] = prob[5]; + x3[4] = x[4]; + xQA3[2] = xQA[2]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[5])/(dedxExp[5]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(x[6] > 0.5){ + xQA[4] = (tof - inttimes[5])/expTOFsigma[5]; // TOF + xQA3[4] = xQA[4]; // TOF + } + if(TMath::Abs(xQA[3]) < 5){ + if(iV0 && fV2) fPhiTracks->Fill(phi,5); + if(fV2) contV0[iV0]->Fill(x,5); + if(fV2) QA[iV0]->Fill(xQA,4); + if(fV3) contV0v3[iV0]->Fill(x3,5); + if(fV3) QAv3[iV0]->Fill(xQA3,4); + } + } + else if(prob[6] > 0.6){ // t + phi[3] = prob[6]; + x[4] = prob[6]; + xQA[2] = prob[6]; + x3[4] = x[4]; + xQA3[2] = xQA[2]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[6])/(dedxExp[6]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(x[6] > 0.5){ + xQA[4] = (tof - inttimes[6])/expTOFsigma[6]; // TOF + xQA3[4] = xQA[4]; // TOF + } + if(TMath::Abs(xQA[3]) < 5){ + if(iV0 && fV2) fPhiTracks->Fill(phi,6); + if(fV2) contV0[iV0]->Fill(x,6); + if(fV2) QA[iV0]->Fill(xQA,5); + if(fV3) contV0v3[iV0]->Fill(x3,6); + if(fV3) QAv3[iV0]->Fill(xQA3,5); + } + } + else if(prob[7] > 0.6){ // He3 + phi[3] = prob[7]; + phi[1] *= 2; + x[4] = prob[7]; + x[3] *= 2; + xQA[2] = prob[7]; + x3[3] = x[3]; + x3[4] = x[4]; + xQA3[2] = xQA[2]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[7])/(dedxExp[7]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(x[6] > 0.5){ + xQA[4] = (tof - inttimes[7])/expTOFsigma[7]; // TOF + xQA3[4] = xQA[4]; // TOF + } + if(TMath::Abs(xQA[3]) < 5){ + if(iV0 && fV2) fPhiTracks->Fill(phi,7); + if(fV2) contV0[iV0]->Fill(x,7); + if(fV2) QA[iV0]->Fill(xQA,6); + if(fV3) contV0v3[iV0]->Fill(x3,7); + if(fV3) QAv3[iV0]->Fill(xQA3,6); + } + x[3] *= 0.5; + x3[3] = x[3]; + } + + if(x[6] > 0.5){ // if TOF was present redo TPC stand alone PID to check the PID in the same acceptance (PID mask = 2) + fPID->ResetDetOR(1); // exclude TOF from PID + tofMismProb = 0; + + fPID->ComputeProb(aodTrack); + probRead = fPID->GetProb(); + + fPID->SetDetOR(1); // include TOF for PID + } + Float_t probTPC[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]}; // TPC stand alone prbabilities + + //pid selection TPC S.A. with TOF matching + x[6]*=2; // set the mask to 2 id TOF is present + if(x[6]<1 || !(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID S.A. PID in TOF acceptance + else if(probTPC[2] > 0.6){ // pi + x[4] = probTPC[2]; + x3[4] = x[4]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[2])/(dedxExp[2]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(TMath::Abs(xQA[3]) < 5){ + if(fV2) contV0[iV0]->Fill(x,1); + if(fV3) contV0v3[iV0]->Fill(x3,1); + } + } + else if(probTPC[3] > 0.6){ // K + x[4] = probTPC[3]; + x3[4] = x[4]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[3])/(dedxExp[3]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(TMath::Abs(xQA[3]) < 5){ + if(fV2) contV0[iV0]->Fill(x,2); + if(fV3) contV0v3[iV0]->Fill(x3,2); + } + } + else if(probTPC[4] > 0.6){ // p + x[4] = probTPC[4]; + x3[4] = x[4]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[4])/(dedxExp[4]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(TMath::Abs(xQA[3]) < 5){ + if(fV2) contV0[iV0]->Fill(x,3); + if(fV3) contV0v3[iV0]->Fill(x3,3); + } + } + else if(probTPC[0] > 0.6){ // e + x[4] = probTPC[0]; + x3[4] = x[4]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[0])/(dedxExp[0]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(TMath::Abs(xQA[3]) < 5){ + if(fV2) contV0[iV0]->Fill(x,4); + if(fV3) contV0v3[iV0]->Fill(x3,4); + } + } + else if(probTPC[5] > 0.6){ // d + x[4] = probTPC[5]; + x3[4] = x[4]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[5])/(dedxExp[5]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(TMath::Abs(xQA[3]) < 5){ + if(fV2) contV0[iV0]->Fill(x,5); + if(fV3) contV0v3[iV0]->Fill(x3,5); + } + } + else if(probTPC[6] > 0.6){ // t + x[4] = probTPC[6]; + x3[4] = x[4]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[6])/(dedxExp[6]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(TMath::Abs(xQA[3]) < 5){ + if(fV2) contV0[iV0]->Fill(x,6); + if(fV3) contV0v3[iV0]->Fill(x3,6); + } + } + else if(probTPC[7] > 0.6){ // He3 + x[4] = probTPC[7]; + x3[4] = x[4]; + x[3] *= 2; + x3[3] = x[3]; + if(dedx > 10.){ + xQA[3] = (dedx - dedxExp[7])/(dedxExp[7]*0.07); // TPC + xQA3[3] = xQA[3]; // TPC + } + if(TMath::Abs(xQA[3]) < 5){ + if(fV2) contV0[iV0]->Fill(x,7); + if(fV3) contV0v3[iV0]->Fill(x3,7); + } + x[3] *= 0.5; + x3[3] = x[3]; + } + } // end side loop + } // end track loop + + // Fill EP distribution histograms + if(fV2) fPhiRPv0A->Fill(iC,evPlAngV0ACor2); + if(fV2) fPhiRPv0C->Fill(iC,evPlAngV0CCor2); + + if(fV3) fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3); + if(fV3) fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3); + + // TPC EP needed for resolution studies (TPC subevent) + Double_t Qx2 = 0, Qy2 = 0; + Double_t Qx3 = 0, Qy3 = 0; + + for(Int_t iT = 0; iT < nAODTracks; iT++) { + + AliAODTrack* aodTrack = aodEvent->GetTrack(iT); + + if (!aodTrack){ + aodTrack->Delete(); + continue; + } + + Bool_t trkFlag = aodTrack->TestFilterBit(1); + + if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || !trkFlag) + continue; + + Double_t b[2] = {-99., -99.}; + Double_t bCov[3] = {-99., -99., -99.}; + if (!aodTrack->PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) + continue; + + if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4)) + continue; + + Qx2 += TMath::Cos(2*aodTrack->Phi()); + Qy2 += TMath::Sin(2*aodTrack->Phi()); + Qx3 += TMath::Cos(3*aodTrack->Phi()); + Qy3 += TMath::Sin(3*aodTrack->Phi()); + + } + + evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.; + evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.; + + // Fill histograms needed for resolution evaluation + if(fV2) fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2))); + if(fV2) fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2))); + if(fV2) fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2))); + + if(fV3) fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3))); + if(fV3) fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3))); + if(fV3) fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3))); + } + +} + +//_____________________________________________________________________________ +Float_t AliAnalysisTaskVnV0::GetVertex(AliAODEvent* aod) const +{ + + Float_t zvtx = -999; + + const AliAODVertex* vtxAOD = aod->GetPrimaryVertex(); + if (!vtxAOD) + return zvtx; + if(vtxAOD->GetNContributors()>0) + zvtx = vtxAOD->GetZ(); + + return zvtx; +} +//_____________________________________________________________________________ +void AliAnalysisTaskVnV0::Terminate(Option_t *) +{ + // Terminate loop + Printf("Terminate()"); +} +//_____________________________________________________________________________ +void AliAnalysisTaskVnV0::OpenInfoCalbration(Int_t run){ + TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root"; + TFile *foadb = TFile::Open(oadbfilename.Data()); + + if(!foadb){ + printf("OADB file %s cannot be opened\n",oadbfilename.Data()); + return; + } + + AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr"); + if(!cont){ + printf("OADB object hMultV0BefCorr is not available in the file\n"); + return; + } + + if(!(cont->GetObject(run))){ + printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run); + run = 137366; + } + fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX(); + + TF1 *fpol0 = new TF1("fpol0","pol0"); + fMultV0->Fit(fpol0,"","",0,31); + fV0Cpol = fpol0->GetParameter(0); + fMultV0->Fit(fpol0,"","",32,64); + fV0Apol = fpol0->GetParameter(0); + + for(Int_t iside=0;iside<2;iside++){ + for(Int_t icoord=0;icoord<2;icoord++){ + for(Int_t i=0;i < nCentrBin;i++){ + char namecont[100]; + if(iside==0 && icoord==0) + sprintf(namecont,"hQxc2_%i",i); + else if(iside==1 && icoord==0) + sprintf(namecont,"hQxa2_%i",i); + else if(iside==0 && icoord==1) + sprintf(namecont,"hQyc2_%i",i); + else if(iside==1 && icoord==1) + sprintf(namecont,"hQya2_%i",i); + + cont = (AliOADBContainer*) foadb->Get(namecont); + if(!cont){ + printf("OADB object %s is not available in the file\n",namecont); + return; + } + + if(!(cont->GetObject(run))){ + printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run); + run = 137366; + } + fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean(); + fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS(); + + //for v3 + if(iside==0 && icoord==0) + sprintf(namecont,"hQxc3_%i",i); + else if(iside==1 && icoord==0) + sprintf(namecont,"hQxa3_%i",i); + else if(iside==0 && icoord==1) + sprintf(namecont,"hQyc3_%i",i); + else if(iside==1 && icoord==1) + sprintf(namecont,"hQya3_%i",i); + + cont = (AliOADBContainer*) foadb->Get(namecont); + if(!cont){ + printf("OADB object %s is not available in the file\n",namecont); + return; + } + + if(!(cont->GetObject(run))){ + printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run); + run = 137366; + } + fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean(); + fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS(); + + } + } + } +} diff --git a/PWG/FLOW/Tasks/AliAnalysisTaskVnV0.h b/PWG/FLOW/Tasks/AliAnalysisTaskVnV0.h new file mode 100644 index 00000000000..80dab1d9253 --- /dev/null +++ b/PWG/FLOW/Tasks/AliAnalysisTaskVnV0.h @@ -0,0 +1,102 @@ +#ifndef ALIANALYSISTASKVNV0_H +#define ALIANALYSISTASKVNV0_H + +// ROOT includes +#include +#include +#include "TTree.h" +#include +#include + +// AliRoot includes +#include +#include +#include +#include "AliFlowBayesianPID.h" + +class TH2F; + +class AliAnalysisTaskVnV0 : public AliAnalysisTaskSE { + public: + AliAnalysisTaskVnV0(); + AliAnalysisTaskVnV0(const char *name); + + virtual ~AliAnalysisTaskVnV0(); + + virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t *option); + virtual void Terminate(Option_t *); + + Double_t GetVtxCut() { return fVtxCut; } + Double_t GetEtaCut() { return fEtaCut; } + Double_t GetMinPt() { return fMinPt; } + + virtual void SetVtxCut(Double_t vtxCut){fVtxCut = vtxCut;} + virtual void SetEtaCut(Double_t etaCut){fEtaCut = etaCut;} + virtual void SetMinPt(Double_t value) {fMinPt = value;} + virtual void SetV2(Bool_t val){fV2 = val;}; + virtual void SetV3(Bool_t val){fV3 = val;}; + + void OpenInfoCalbration(Int_t run); + + private: + AliAnalysisTaskVnV0(const AliAnalysisTaskVnV0 &); + AliAnalysisTaskVnV0 & operator=(const AliAnalysisTaskVnV0 &); + + virtual Float_t GetVertex(AliAODEvent* aod) const; + virtual void Analyze(AliAODEvent* aodEvent, Float_t v0Centr); + + AliAODEvent* fAOD; //! AOD object + + static const Int_t nCentrBin = 9; // # cenrality bins + + // + // Cuts and options + // + Double_t fVtxCut; // Vtx cut on z position in cm + Double_t fEtaCut; // Eta cut used to select particles + Double_t fMinPt; // Min pt - for histogram limits + + Int_t fRun; // current run checked to load VZERO calibrations + + TList *fList,*fList2; // List for output objects + // + // Output objects + TProfile *fMultV0; // object containing VZERO calibration information + Float_t fV0Cpol,fV0Apol; // loaded by OADB + Float_t fMeanQ[nCentrBin][2][2]; // and recentering + Float_t fWidthQ[nCentrBin][2][2]; // ... + Float_t fMeanQv3[nCentrBin][2][2]; // also for v3 + Float_t fWidthQv3[nCentrBin][2][2]; // ... + + AliCFContainer *fContAllChargesV0A; // containers for v2 (A nd C) and v3 (A and C) + AliCFContainer *fContAllChargesV0C; // ... + AliCFContainer *fContAllChargesV0Av3; // ... + AliCFContainer *fContAllChargesV0Cv3; // ... + + TProfile *fHResTPCv0A2,*fHResTPCv0C2,*fHResv0Cv0A2; // TProfile for subevent resolution (output) + TProfile *fHResTPCv0A3,*fHResTPCv0C3,*fHResv0Cv0A3; // also for v3 + + TH2F *fPhiRPv0A,*fPhiRPv0C; // EP distribution vs. centrality (v2) + TH2F *fPhiRPv0Av3,*fPhiRPv0Cv3; // EP distribution vs. centrality (v3) + + AliCFContainer *fPhiTracks; // phi distribution of particles (if needed) + + + AliCFContainer *fQA,*fQA2; // QA container (v2) + AliCFContainer *fQAv3,*fQA2v3; // QA container (v3) + + AliFlowBayesianPID *fPID; // PID class for the Bayesian probabilities + + TTree *fTree; // tree to debug EP (if needed) + + Float_t fCentrality; // current centrality for the tree + Float_t evPlAngV0ACor2,evPlAngV0CCor2,evPlAng2; // subevent EPs (v2) + Float_t evPlAngV0ACor3,evPlAngV0CCor3,evPlAng3; // subevent EPs (v3) + + Bool_t fV2,fV3; // swith to set the armonics + + ClassDef(AliAnalysisTaskVnV0, 1); //Analysis task v2 and v3 analysis on AOD +}; + +#endif diff --git a/PWG/PWGflowTasksLinkDef.h b/PWG/PWGflowTasksLinkDef.h index f342bcde0e8..d6648e663c8 100644 --- a/PWG/PWGflowTasksLinkDef.h +++ b/PWG/PWGflowTasksLinkDef.h @@ -29,6 +29,7 @@ #pragma link C++ class AliFlowBayesianPID+; #pragma link C++ class AliAnalysisTaskPhiFlow+; #pragma link C++ class AliAnalysisTaskFilterFE+; +#pragma link C++ class AliAnalysisTaskVnV0+; #endif diff --git a/PWGCF/FLOW/macros/AddTaskVZERO.C b/PWGCF/FLOW/macros/AddTaskVZERO.C new file mode 100644 index 00000000000..43b1180e91a --- /dev/null +++ b/PWGCF/FLOW/macros/AddTaskVZERO.C @@ -0,0 +1,45 @@ +AliAnalysisTask *AddTaskVZERO(AliAnalysisManager *mgr,Bool_t kV2=kTRUE,Bool_t kV3=kTRUE){ + char fileout[100]; + sprintf(fileout,"outVZEROv2.root"); + char fileout2[100]; + sprintf(fileout2,"outVZEROv3.root"); + + //get the current analysis manager + // AnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error("No manager found in AddTaskVZERO. Why?"); + return 0; + } + // currently don't accept AOD input + if (!mgr->GetInputEventHandler()->InheritsFrom(AliAODInputHandler::Class())) { + Error("AddTaskVZERO","This task works only with AOD input!"); + return 0; + } + + //========= Add tender to the ANALYSIS manager and set default storage ===== + char mytaskName[100]; + sprintf(mytaskName,"AliAnalysisTaskVnV0.cxx"); + + AliAnalysisTaskVnV0 *task = new AliAnalysisTaskVnV0(mytaskName); + task->SetV2(kV2); + task->SetV3(kV3); + + mgr->AddTask(task); + + //Attach input to my tasks + AliAnalysisDataContainer *cinput = mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer); + mgr->ConnectInput(task,0,mgr->GetCommonInputContainer()); + + // Attach output to my tasks + if(kV2){ + AliAnalysisDataContainer *cOutputL= mgr->CreateContainer("contVZEROv2",TList::Class(), AliAnalysisManager::kOutputContainer, fileout); + mgr->ConnectOutput(task, 1, cOutputL); + } + if(kV3){ + AliAnalysisDataContainer *cOutputL2= mgr->CreateContainer("contVZEROv3",TList::Class(), AliAnalysisManager::kOutputContainer, fileout2); + mgr->ConnectOutput(task, 2, cOutputL2); + } + printf("task really added\n"); + + return task; +} diff --git a/PWGCF/FLOW/macros/extractFlowVZERO.C b/PWGCF/FLOW/macros/extractFlowVZERO.C new file mode 100644 index 00000000000..94ec7c179df --- /dev/null +++ b/PWGCF/FLOW/macros/extractFlowVZERO.C @@ -0,0 +1,279 @@ +extractFlowVZERO(Int_t itech=0,Int_t ibin = 3,Int_t step = 0,Int_t arm=2,Float_t pTh = 0.9,Int_t charge=0,Int_t addbins=0){ + TCanvas *can1 = new TCanvas("cV0Acheck","cV0Acheck"); + +/* +read outVZEROv2.root or outVZEROv3.root files + +itech=0 TPC stand alone PID and selection +itech=1 TPC stand alone PID but TOF track selection +itech=2 TPC&TOF PID (TOF strictly required) +itech=3 TPC|TOF PID (the same of before but using TPC stand alone when TOF not available) + +ibin = centrality bin + +step = species (0=all charges, 1=pi, 2=K, 3=p, 4=e, 5=d, 6=t, 7=3He + +arm = armonic 2 (v2) or 3 (v3) + +pTh = probability threshold = 0.6 , 0.8 , 0.9 , 0.95 + +charge = 1(pos), -1(neg), 0(both) + +addbins = to merge more centrality bin (i.e. ibin = 2, addbins = 3 merge 2,3,4,5 centrality bins) - preliminary version!!!!!! + +*/ + + // load libraries + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libCORRFW.so"); + + // get objects + char filein[100]; + sprintf(filein,"outVZEROv%i.root",arm); + TFile *f=new TFile(filein); + sprintf(filein,"contVZEROv%i",arm); + TList *l = f->Get(filein); + l->ls(); + AliCFContainer *c = l->At(0); + AliCFContainer *c3 = l->At(1); + TProfile *p1 = l->At(2); + TProfile *p2 = l->At(3); + TProfile *p3 = l->At(4); + Int_t iVar[2] = {1,3}; + + Int_t iVarPsi[1] = {5}; + /* 0=centrality , 1= cos(arm*deltaphi) , 2=charge , 3=pt , 4=probability , 5=EP , 6=PIDmask*/ + Double_t iMin[7] = {ibin,-1,-1.5,0,pTh+0.00001,-3.14,0}; + Double_t iMax[7] = {ibin+addbins,1,1.5,20,1.1,3.14,1}; + + char *iTechName[4] = {"tpc","tpcInTof","tof","tpcORtof"}; + + Bool_t kOnlyTPC = kFALSE; + + if(itech==0 && (step!=0)){ + kOnlyTPC = kTRUE; + } + else if(itech==1 && (step!=0)){ + iMin[6] = 2; + iMax[6] = 2; + } + else if(itech==2 || itech==1){ + iMin[6] = 1; + iMax[6] = 1; + } + else{ + iMin[6] = 0; + iMax[6] = 1; + } + + if(charge==-1) iMax[2]=0; + if(charge==1) iMin[2]=0; + + // EP distribution needed for NUA corrections + TH2F *hPhiA = l->At(5); + TH2F *hPhiC = l->At(6); + + TH1D *hPhiPA = hPhiA->ProjectionY("_py",iMin[0]+1,iMax[0]+1); + TH1D *hPhiPC = hPhiC->ProjectionY("_py",iMin[0]+1,iMax[0]+1); + + // EP resolution variables + Float_t res1=0,res2=0,res3=0; + Float_t eres1=0,eres2=0,eres3=0; + + for(Int_t i=iMin[0];i<=iMax[0];i++){ // in case of more centrality bins weighting with the errors + if(p1->GetBinError(i+1)){ + eres1 += 1./p1->GetBinError(i+1)/p1->GetBinError(i+1); + res1 += p1->GetBinContent(i+1)/p1->GetBinError(i+1)/p1->GetBinError(i+1); + } + if(p2->GetBinError(i+1)){ + eres2 += 1./p2->GetBinError(i+1)/p2->GetBinError(i+1); + res2 += p2->GetBinContent(i+1)/p2->GetBinError(i+1)/p2->GetBinError(i+1); + } + if(p3->GetBinError(i+1)){ + eres3 += 1./p3->GetBinError(i+1)/p3->GetBinError(i+1); + res3 += p3->GetBinContent(i+1)/p3->GetBinError(i+1)/p3->GetBinError(i+1); + } + + if(eres1) res1 /= eres1; + if(eres2) res2 /= eres2; + if(eres3) res3 /= eres3; + + if(eres1) eres1 = TMath::Sqrt(1./eres1); + if(eres1) eres2 = TMath::Sqrt(1./eres2); + if(eres1) eres3 = TMath::Sqrt(1./eres3); + } + + // NUA correction (fit to EP distribution) + TF1 *fpol = new TF1("fPol","pol0"); + hPhiPA->Fit("fPol","","",-TMath::Pi()/arm,TMath::Pi()/arm); + Float_t scalA = fPol->GetParameter(0); + hPhiPC->Fit("fPol","","",-TMath::Pi()/arm,TMath::Pi()/arm); + Float_t scalC = fPol->GetParameter(0); + + AliCFContainer *c2[20]; + AliCFContainer *c4[20]; + TH2D *h[20],*h2[20]; + + AliCFContainer *c2bis[20]; + AliCFContainer *c4bis[20]; + + AliCFContainer *cPsi2[20]; + AliCFContainer *cPsi4[20]; + TH1D *hPsi[20],*hPsi2[20]; + AliCFContainer *cPsi2bis[20]; + AliCFContainer *cPsi4bis[20]; + + Float_t intA = 0; + Float_t intC = 0; + + for(Int_t i=5;i<15;i++){ + printf("%i\n",i); + iMin[5] = hPhiPA->GetBinCenter(i+1); + iMax[5] = hPhiPA->GetBinCenter(i+1); + if(!kOnlyTPC){ + c2[i] = c->MakeSlice(2,iVar,iMin,iMax); + c4[i]= c3->MakeSlice(2,iVar,iMin,iMax); + } + else{ // merge to maskPID bins (TPC standalone without TOF + TPC stand alone with TOF) + iMin[6] = 0; + iMax[6] = 0; + c2[i] = c->MakeSlice(2,iVar,iMin,iMax); + c4[i]= c3->MakeSlice(2,iVar,iMin,iMax); + + iMin[6] = 2; + iMax[6] = 2; + c2bis[i] = c->MakeSlice(2,iVar,iMin,iMax); + c4bis[i]= c3->MakeSlice(2,iVar,iMin,iMax); + } + + + if(!kOnlyTPC){ + cPsi2[i] = c->MakeSlice(1,iVarPsi,iMin,iMax); + cPsi4[i]= c3->MakeSlice(1,iVarPsi,iMin,iMax); + } + else{ // merge to maskPID bins (TPC standalone without TOF + TPC stand alone with TOF) + iMin[6] = 0; + iMax[6] = 0; + cPsi2[i] = c->MakeSlice(1,iVarPsi,iMin,iMax); + cPsi4[i]= c3->MakeSlice(1,iVarPsi,iMin,iMax); + + iMin[6] = 2; + iMax[6] = 2; + cPsi2bis[i] = c->MakeSlice(1,iVarPsi,iMin,iMax); + cPsi4bis[i]= c3->MakeSlice(1,iVarPsi,iMin,iMax); + } + + h[i] = (TH2D *) c2[i]->Project(step,0,1); + h2[i] = (TH2D *) c4[i]->Project(step,0,1); + if(kOnlyTPC){ + TH2D *htemp = (TH2D *) c2bis[i]->Project(step,0,1); + h[i]->Add(htemp); + htemp = (TH2D *) c4bis[i]->Project(step,0,1); + h2[i]->Add(htemp); + } + + if(hPhiPA->GetBinContent(i)) h[i]->Scale(scalA/hPhiPA->GetBinContent(i+1)); // reweighting for NUA correction + if(hPhiPC->GetBinContent(i)) h2[i]->Scale(scalC/hPhiPC->GetBinContent(i+1)); // reweighting for NUA correction + + hPsi[i] = (TH1D *) cPsi2[i]->Project(step,0); + hPsi2[i] = (TH1D *) cPsi4[i]->Project(step,0); + + if(kOnlyTPC){ // merge to maskPID bins (TPC standalone without TOF + TPC stand alone with TOF) + TH1D *htemp2 = (TH1D *) cPsi2bis[i]->Project(step,0); + hPsi[i]->Add(htemp2); + htemp2 = (TH1D *) cPsi4bis[i]->Project(step,0); + hPsi2[i]->Add(htemp2); + } + + if(hPhiPA->GetBinContent(i+1)) hPsi[i]->Scale(scalA/hPhiPA->GetBinContent(i+1)); // check of reweighting for NUA correction + if(hPhiPC->GetBinContent(i+1)) hPsi2[i]->Scale(scalC/hPhiPC->GetBinContent(i+1)); // check of reweighting for NUA correction + + intA += hPsi[i]->Integral(); + intC += hPsi2[i]->Integral(); + + } + +// hPhiPA->Scale(intA / (scalA)/10); +// hPhiPC->Scale(intC / (scalC)/10); + + // NUA correction check V0A + hPhiPA->Draw(); + hPsi[5]->Draw("SAME"); + hPsi[5]->Scale(scalA / intA * 10); + for(Int_t i=6;i<15;i++){ + h[5]->Add(h[i]); + h2[5]->Add(h2[i]); + hPsi[i]->Draw("SAME"); + hPsi[i]->Scale(scalA / intA * 10); + } + + // NUA correction check V0C + TCanvas *can2 = new TCanvas("cV0Ccheck","cV0Ccheck"); + hPhiPC->Draw(); + hPsi2[5]->Draw("SAME"); + for(Int_t i=6;i<15;i++){ + h[5]->Add(h[i]); + h2[5]->Add(h2[i]); + hPsi2[i]->Draw("SAME"); + hPsi2[i]->Scale(scalC /intC * 10); + } + hPsi2[5]->Scale(scalC /intC * 10); + + // Flow for V0A and V0C separately + TCanvas *can3 = new TCanvas("cFlowComp","cFlowComp"); + TProfile *pp = h[5]->ProfileY(); + pp->Draw(); + TProfile *pp2 = h2[5]->ProfileY(); + pp2->Draw("SAME"); + + printf("nev (selected within 0-80%) = %i\n",p1->GetEntries()); + + + // correction for resoltion + Float_t scaling = sqrt(res1*res3/res2); + pp->Scale(1./scaling); + printf("resolution V0A = %f\n",scaling); + Float_t err1_2 = eres1*eres1/res1/res1/4 + + eres2*eres2/res2/res2/4 + + eres3*eres3/res3/res3/4; + Float_t err2_2 = err1_2; + err1_2 /= scaling*scaling; + scaling = sqrt(res2*res3/res1); + err2_2 /= scaling*scaling; + pp2->Scale(1./scaling); + printf("resolution V0C =%f\n",scaling); + + char title[100]; + sprintf(title,"VZERO EP;p_{t} (GeV/c);v_{%i}",arm); + pp->SetTitle(title); + + // Average V0A-V0C + TH1D *pAll = pp->ProjectionX(); + + for(Int_t i=1;i <= pAll->GetNbinsX();i++){ + Float_t e1 = err1_2*pp->GetBinContent(i)*pp->GetBinContent(i) + pp->GetBinError(i)*pp->GetBinError(i); + Float_t e2 = err2_2*pp2->GetBinContent(i)*pp2->GetBinContent(i) + pp2->GetBinError(i)*pp2->GetBinError(i); + Float_t xval = 0,exval = 0; + if(e1 >0 && e2>0){ + xval = (pp->GetBinContent(i)/e1 + pp2->GetBinContent(i)/e2)/(1/e1 + 1/e2); + exval = 1./sqrt(1/e1 + 1/e2); + } + pAll->SetBinContent(i,xval); + pAll->SetBinError(i,exval); + } + // combined measurement + TCanvas *can4 = new TCanvas("cFlowCombined","cFlowCombined"); + pAll->Draw(); + + char name[100]; + sprintf(name,"out%i-%i_%i_%4.2f_%i%sv%i.root",iMin[0],iMax[0],step,pTh,charge,iTechName[itech],arm); + TFile *fout = new TFile(name,"RECREATE"); + + pAll->SetName("histo"); + pAll->Write(); + can1->Write(); + can2->Write(); + can3->Write(); + fout->Close(); +} -- 2.39.3