fContAllChargesMCA(NULL),
fContAllChargesMCC(NULL),
fContAllChargesMCAv3(NULL),
- fContAllChargesMCCv3(NULL)
+ fContAllChargesMCCv3(NULL),
+ fFillDCA(kFALSE),
+ fContQApid(NULL),
+ fModulationDEDx(kFALSE)
{
// Default constructor (should not be used)
fList->SetName("resultsV2");
fContAllChargesMCA(NULL),
fContAllChargesMCC(NULL),
fContAllChargesMCAv3(NULL),
- fContAllChargesMCCv3(NULL)
+ fContAllChargesMCCv3(NULL),
+ fFillDCA(kFALSE),
+ fContQApid(NULL),
+ fModulationDEDx(kFALSE)
{
DefineOutput(1, TList::Class());
const Int_t nPsiTOFres = 10;
const Int_t nMaskPID = 3;
- Int_t binsTOF[5] = {nCentrTOFres,nChargeBinsTOFres,nProbTOFres,nPsiTOFres,nMaskPID};
+ Int_t nDCABin = 1; // put to 1 not to store this info
+ if(fFillDCA) nDCABin = 3;
+ if(fIsMC && nDCABin>1) nDCABin = 6;
+ /*
+ 0 = DCAxy < 2.4 && all (or Physical primary if MC)
+ 1 = DCAxy > 2.4 && all (or Physical primary if MC)
+ 2 = DCAxy < 2.4 && not Physical Primary for MC
+ 3 = DCAxy > 2.4 && not Physical Primary for MC
+ */
+
+ Int_t binsTOF[6] = {nCentrTOFres,nChargeBinsTOFres,nProbTOFres,nPsiTOFres,nMaskPID,nDCABin};
Int_t binsTOFmc[5] = {nCentrTOFres,nChargeBinsTOFres,1,nPsiTOFres,2};
Int_t binsTOFmcPureMC[5] = {nCentrTOFres,nChargeBinsTOFres,1,nPsiTOFres,1};
// v2 container
- fContAllChargesV0A = new AliFlowVZEROResults("v2A",5,binsTOF);
+ fContAllChargesV0A = new AliFlowVZEROResults("v2A",6,binsTOF);
fContAllChargesV0A->SetVarRange(0,-0.5,8.5); // centrality
fContAllChargesV0A->SetVarRange(1,-1.5,1.5); // charge
fContAllChargesV0A->SetVarRange(2,0.6,1.0001);// prob
fContAllChargesV0A->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
fContAllChargesV0A->SetVarRange(4,-0.5,2.5); // pid mask
+ fContAllChargesV0A->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
fContAllChargesV0A->SetVarName(0,"centrality");
fContAllChargesV0A->SetVarName(1,"charge");
fContAllChargesV0A->SetVarName(2,"prob");
fContAllChargesV0A->SetVarName(3,"#Psi");
fContAllChargesV0A->SetVarName(4,"PIDmask");
+ fContAllChargesV0A->SetVarName(5,"DCAbin");
if(fV2) fContAllChargesV0A->AddSpecies("all",nPtBinsTOF,binsPtTOF);
if(fV2) fContAllChargesV0A->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
if(fV2) fContAllChargesV0A->AddSpecies("k",nPtBinsTOF,binsPtTOF);
if(fV2) fContAllChargesV0A->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
if(fV2) fContAllChargesV0A->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
- fContAllChargesV0C = new AliFlowVZEROResults("v2C",5,binsTOF);
+ fContAllChargesV0C = new AliFlowVZEROResults("v2C",6,binsTOF);
fContAllChargesV0C->SetVarRange(0,-0.5,8.5); // centrality
fContAllChargesV0C->SetVarRange(1,-1.5,1.5); // charge
fContAllChargesV0C->SetVarRange(2,0.6,1.0001);// prob
fContAllChargesV0C->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
fContAllChargesV0C->SetVarRange(4,-0.5,2.5); // pid mask
+ fContAllChargesV0C->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
fContAllChargesV0C->SetVarName(0,"centrality");
fContAllChargesV0C->SetVarName(1,"charge");
fContAllChargesV0C->SetVarName(2,"prob");
fContAllChargesV0C->SetVarName(3,"#Psi");
fContAllChargesV0C->SetVarName(4,"PIDmask");
+ fContAllChargesV0C->SetVarName(5,"DCAbin");
if(fV2) fContAllChargesV0C->AddSpecies("all",nPtBinsTOF,binsPtTOF);
if(fV2) fContAllChargesV0C->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
if(fV2) fContAllChargesV0C->AddSpecies("k",nPtBinsTOF,binsPtTOF);
}
// v3 container
- fContAllChargesV0Av3 = new AliFlowVZEROResults("v3A",5,binsTOF);
+ fContAllChargesV0Av3 = new AliFlowVZEROResults("v3A",6,binsTOF);
fContAllChargesV0Av3->SetVarRange(0,-0.5,8.5); // centrality
fContAllChargesV0Av3->SetVarRange(1,-1.5,1.5); // charge
fContAllChargesV0Av3->SetVarRange(2,0.6,1.0001);// prob
fContAllChargesV0Av3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
+ fContAllChargesV0Av3->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
fContAllChargesV0Av3->SetVarRange(4,-0.5,2.5); // pid mask
fContAllChargesV0Av3->SetVarName(0,"centrality");
fContAllChargesV0Av3->SetVarName(1,"charge");
fContAllChargesV0Av3->SetVarName(2,"prob");
fContAllChargesV0Av3->SetVarName(3,"#Psi");
fContAllChargesV0Av3->SetVarName(4,"PIDmask");
+ fContAllChargesV0Av3->SetVarName(5,"DCAbin");
if(fV3) fContAllChargesV0Av3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
if(fV3) fContAllChargesV0Av3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
if(fV3) fContAllChargesV0Av3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
if(fV3) fContAllChargesV0Av3->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
if(fV3) fContAllChargesV0Av3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
- fContAllChargesV0Cv3 = new AliFlowVZEROResults("v3C",5,binsTOF);
+ fContAllChargesV0Cv3 = new AliFlowVZEROResults("v3C",6,binsTOF);
fContAllChargesV0Cv3->SetVarRange(0,-0.5,8.5); // centrality
fContAllChargesV0Cv3->SetVarRange(1,-1.5,1.5); // charge
fContAllChargesV0Cv3->SetVarRange(2,0.6,1.0001);// prob
fContAllChargesV0Cv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
fContAllChargesV0Cv3->SetVarRange(4,-0.5,2.5); // pid mask
+ fContAllChargesV0Cv3->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
fContAllChargesV0Cv3->SetVarName(0,"centrality");
fContAllChargesV0Cv3->SetVarName(1,"charge");
fContAllChargesV0Cv3->SetVarName(2,"prob");
fContAllChargesV0Cv3->SetVarName(3,"#Psi");
fContAllChargesV0Cv3->SetVarName(4,"PIDmask");
+ fContAllChargesV0Cv3->SetVarName(5,"DCAbin");
if(fV3) fContAllChargesV0Cv3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
if(fV3) fContAllChargesV0Cv3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
if(fV3) fContAllChargesV0Cv3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
// fList->Add(fTree); // comment if not needed
+ const Int_t nDCA = 300;
+ Double_t DCAbin[nDCA+1];
+ for(Int_t i=0;i <= nDCA;i++){
+ DCAbin[i] = -3 +i*6.0/nDCA;
+ }
+
+ char nameHistos[100];
+ for(Int_t iC=0;iC < nCentrBin;iC++){
+ snprintf(nameHistos,100,"fHdcaPtPiCent%i",iC);
+ fHdcaPt[iC][0] = new TH2D(nameHistos,"DCA_{xy} for #pi;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtKaCent%i",iC);
+ fHdcaPt[iC][1] = new TH2D(nameHistos,"DCA_{xy} for K;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtPrCent%i",iC);
+ fHdcaPt[iC][2] = new TH2D(nameHistos,"DCA_{xy} for #bar{p};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtElCent%i",iC);
+ fHdcaPt[iC][3] = new TH2D(nameHistos,"DCA_{xy} for e;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtDeCent%i",iC);
+ fHdcaPt[iC][4] = new TH2D(nameHistos,"DCA_{xy} for #bar{d};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtTrCent%i",iC);
+ fHdcaPt[iC][5] = new TH2D(nameHistos,"DCA_{xy} for #bar{t};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtHeCent%i",iC);
+ fHdcaPt[iC][6] = new TH2D(nameHistos,"DCA_{xy} for #bar{^{3}He};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ }
+
+ if(fFillDCA && fQAsw){
+ for(Int_t i=0;i<7;i++)
+ for(Int_t iC=0;iC < nCentrBin;iC++)
+ fList4->Add(fHdcaPt[iC][i]);
+ }
+ if(fIsMC){
+ for(Int_t iC=0;iC < nCentrBin;iC++){
+ snprintf(nameHistos,100,"fHdcaPtPiSecCent%i",iC);
+ fHdcaPtSec[iC][0] = new TH2D(nameHistos,"DCA_{xy} for secondary #pi;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtKaSecCent%i",iC);
+ fHdcaPtSec[iC][1] = new TH2D(nameHistos,"DCA_{xy} for secondary K;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtPrSecCent%i",iC);
+ fHdcaPtSec[iC][2] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{p};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtElSecCent%i",iC);
+ fHdcaPtSec[iC][3] = new TH2D(nameHistos,"DCA_{xy} for secondary e;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtDeSecCent%i",iC);
+ fHdcaPtSec[iC][4] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{d};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtTrSecCent%i",iC);
+ fHdcaPtSec[iC][5] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{t};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ snprintf(nameHistos,100,"fHdcaPtHeSecCent%i",iC);
+ fHdcaPtSec[iC][6] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{^{3}He};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
+ }
+
+ if(fFillDCA && fQAsw){
+ for(Int_t i=0;i<7;i++)
+ for(Int_t iC=0;iC < nCentrBin;iC++)
+ fList4->Add(fHdcaPtSec[iC][i]);
+ }
+ }
+
+ // Add TProfile Extra QA
+ const Int_t nBinQApid = 2;
+ Int_t binQApid[nBinQApid] = {nCentrTOF,200};
+ const Int_t nbinsigma = 100;
+ Double_t nsigmaQA[nbinsigma+1];
+ for(Int_t i=0;i<nbinsigma+1;i++){
+ nsigmaQA[i] = -10 + 20.0*i/nbinsigma;
+ }
+ fContQApid = new AliFlowVZEROResults("qaPID",nBinQApid,binQApid);
+ fContQApid->SetVarRange(0,-0.5,8.5); // centrality
+ fContQApid->SetVarRange(1,0,20); // charge
+ fContQApid->SetVarName(0,"centrality");
+ fContQApid->SetVarName(1,"p_{t}");
+ fContQApid->AddSpecies("piTPC",nbinsigma,nsigmaQA);
+ fContQApid->AddSpecies("piTOF",nbinsigma,nsigmaQA);
+ fContQApid->AddSpecies("kaTPC",nbinsigma,nsigmaQA);
+ fContQApid->AddSpecies("kaTOF",nbinsigma,nsigmaQA);
+ fContQApid->AddSpecies("prTPC",nbinsigma,nsigmaQA);
+ fContQApid->AddSpecies("prTOF",nbinsigma,nsigmaQA);
+ if(fV2) fList->Add(fContQApid);
+
printf("Output creation ok!!\n\n\n\n");
// Post output data.
}
Bool_t trkFlag = aodTrack->TestFilterBit(1); // TPC only tracks
+ if(fFillDCA) trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < 70) || !trkFlag){
continue;
if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
continue;
- if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
+ if (!fFillDCA && ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4)))
continue;
+ if(fFillDCA && TMath::Abs(b[0]) > 3.0 && TMath::Abs(b[1]) > 3)
+ 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};
AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
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)
+ 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)
Float_t v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
Float_t v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
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
- Float_t x[5] = {iC,aodTrack->Charge(),1,evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5}; // to fill analysis v2 container
- Float_t x3[5] = {iC,aodTrack->Charge(),1,evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5}; // to fill analysis v3 container
+ Float_t x[6] = {iC,aodTrack->Charge(),1,evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v2 container
+ Float_t x3[6] = {iC,aodTrack->Charge(),1,evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v3 container
+
+ // in case fill DCA info
+ if(fFillDCA){
+ if(TMath::Abs(b[0]) > 0.1){
+ x[5] = 1;
+ x3[5] = 1;
+ }
+ if(TMath::Abs(b[0]) > 0.3){
+ x[5] = 2;
+ x3[5] = 2;
+ }
+ if(fIsMC && mcArray){
+ if(!((AliAODMCParticle*)mcArray->At(TMath::Abs(aodTrack->GetLabel())))->IsPhysicalPrimary()){
+ x[5] += 3;
+ x3[5] += 3;
+ }
+ }
+ }
// Fill no PID
if(fV2) contV0[iV0]->Fill(0,aodTrack->Pt(),v2V0,x);
Float_t xQA[5] = {iC,aodTrack->Pt(), 0.0,deltaPhiV0,x[4]}; // v2
Float_t xQA3[5] = {iC,aodTrack->Pt(), 0.0,deltaPhiV0v3,x[4]}; // v3
+ // extra QA TProfiles
+ if(iV0==1 && aodTrack->Pt() < 20 && fPID->GetCurrentMask(0) && fPID->GetCurrentMask(1)){
+ Float_t xQApid[2] = {iC,aodTrack->Pt()};
+ fContQApid->Fill(0,nsigmaTPC[2],v2V0,xQApid); // v2 TPC (V0C) w.r.t pions
+ fContQApid->Fill(1,nsigmaTOF[2],v2V0,xQApid); // v2 TOF (V0C) w.r.t. pions
+ fContQApid->Fill(2,nsigmaTPC[3],v2V0,xQApid); // v2 TPC (V0C) w.r.t kaons
+ fContQApid->Fill(3,nsigmaTOF[3],v2V0,xQApid); // v2 TOF (V0C) w.r.t. kaons
+ fContQApid->Fill(4,nsigmaTPC[4],v2V0,xQApid); // v2 TPC (V0C) w.r.t protons
+ fContQApid->Fill(5,nsigmaTOF[4],v2V0,xQApid); // v2 TOF (V0C) w.r.t. protons
+ }
+
// QA fill
if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid() || dedx < 10. || aodTrack->Pt() < 0 || aodTrack->Pt() > 7){}
else{
if(TMath::Abs(nsigmaTPC[2]) < 5){ // TPC 5 sigma extra cut to accept the track
if(fV2) contV0[iV0]->Fill(1,aodTrack->Pt(),v2V0,x);
if(fV3) contV0v3[iV0]->Fill(1,aodTrack->Pt(),v3V0,x3);
+ if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][0]->Fill(aodTrack->Pt(),b[0]);
+ else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][0]->Fill(aodTrack->Pt(),b[0]);
}
}
else if(prob[3] > 0.6){ // K
if(TMath::Abs(nsigmaTPC[3]) < 5){
if(fV2) contV0[iV0]->Fill(2,aodTrack->Pt(),v2V0,x);
if(fV3) contV0v3[iV0]->Fill(2,aodTrack->Pt(),v3V0,x3);
+ if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][1]->Fill(aodTrack->Pt(),b[0]);
+ else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][1]->Fill(aodTrack->Pt(),b[0]);
}
}
else if(prob[4] > 0.6){ // p
if(TMath::Abs(nsigmaTPC[4]) < 5){
if(fV2) contV0[iV0]->Fill(3,aodTrack->Pt(),v2V0,x);
if(fV3) contV0v3[iV0]->Fill(3,aodTrack->Pt(),v3V0,x3);
+ if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][2]->Fill(aodTrack->Pt(),b[0]);
+ else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][2]->Fill(aodTrack->Pt(),b[0]);
}
}
else if(prob[0] > 0.6){ // e
if(TMath::Abs(nsigmaTPC[0]) < 5){
if(fV2) contV0[iV0]->Fill(4,aodTrack->Pt(),v2V0,x);
if(fV3) contV0v3[iV0]->Fill(4,aodTrack->Pt(),v3V0,x3);
+ if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][3]->Fill(aodTrack->Pt(),b[0]);
+ else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][3]->Fill(aodTrack->Pt(),b[0]);
}
}
else if(prob[1] > 0.6){ // mu
if(TMath::Abs(nsigmaTPC[5]) < 5){
if(fV2) contV0[iV0]->Fill(5,aodTrack->Pt(),v2V0,x);
if(fV3) contV0v3[iV0]->Fill(5,aodTrack->Pt(),v3V0,x3);
+ if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][4]->Fill(aodTrack->Pt(),b[0]);
+ else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][4]->Fill(aodTrack->Pt(),b[0]);
}
}
else if(prob[6] > 0.6){ // t
if(TMath::Abs(nsigmaTPC[6]) < 5){
if(fV2) contV0[iV0]->Fill(6,aodTrack->Pt(),v2V0,x);
if(fV3) contV0v3[iV0]->Fill(6,aodTrack->Pt(),v3V0,x3);
+ if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][5]->Fill(aodTrack->Pt(),b[0]);
+ else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][5]->Fill(aodTrack->Pt(),b[0]);
}
}
else if(prob[7] > 0.6){ // He3
if(TMath::Abs(nsigmaTPC[7]) < 5){
if(fV2) contV0[iV0]->Fill(7,aodTrack->Pt()*2,v2V0,x);
if(fV3) contV0v3[iV0]->Fill(7,aodTrack->Pt()*2,v3V0,x3);
+ if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][6]->Fill(aodTrack->Pt(),b[0]);
+ else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][6]->Fill(aodTrack->Pt(),b[0]);
}
}
void OpenInfoCalbration(Int_t run);
+ void SetFillDCAinfo(Bool_t flag=kTRUE){fFillDCA = flag;};
+
+ void SetModulationDEDx(Bool_t flag=kTRUE){fModulationDEDx=flag;};
+
private:
AliAnalysisTaskVnV0(const AliAnalysisTaskVnV0 &old);
AliAnalysisTaskVnV0& operator=(const AliAnalysisTaskVnV0 &source);
AliFlowVZEROResults *fContAllChargesMCAv3; //! results
AliFlowVZEROResults *fContAllChargesMCCv3; //! results
- ClassDef(AliAnalysisTaskVnV0, 5); //Analysis task v2 and v3 analysis on AOD
+ Bool_t fFillDCA; // require to fill also DCA info
+ TH2D *fHdcaPt[nCentrBin][7]; //! DCA distribution (for MC primary)
+ TH2D *fHdcaPtSec[nCentrBin][7]; //! DCA distribution (for MC secondary, not used for data)
+ AliFlowVZEROResults *fContQApid; //! QA pid object
+
+ Bool_t fModulationDEDx; //add a modulation on the dE/dx response w.r.t. EP (kFALSE default)
+
+ ClassDef(AliAnalysisTaskVnV0, 6); //Analysis task v2 and v3 analysis on AOD
};
#endif
ClassImp(AliFlowBayesianPID)
-TH2D* AliFlowBayesianPID::fghPriors[fgkNspecies] = {NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL}; // histo with priors (hardcoded)
+TH2D* AliFlowBayesianPID::fghPriors[fgkNspecies] = {NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL}; // histo with priors (hardcoded)
TSpline3* AliFlowBayesianPID::fgMism = NULL; // function for mismatch
TH1D* AliFlowBayesianPID::fgHtofChannelDist=NULL;
fghPriors[7] = new TH2D("fghPriorsHe","Priors as a function of Centrality [p];centrality;p_{t} (GeV/c)",12,-10,110,80,0,20);
redopriors = kTRUE;
}
+ if(! fghPriors[8]){
+ fghPriors[8] = new TH2D("fghPriorsAlfa","Priors as a function of Centrality [p];centrality;p_{t} (GeV/c)",12,-10,110,80,0,20);
+ redopriors = kTRUE;
+ }
if(redopriors) SetPriors();
fProb[5]=0.0;
fProb[6]=0.0;
fProb[7]=0.0;
+ fProb[8]=0.0;
fMass[0] = fDB->GetParticle(11)->Mass(); // e mass
fMass[1] = fDB->GetParticle(13)->Mass(); // mu mass
fMass[5] = fDB->GetParticle(2212)->Mass()+fDB->GetParticle(2112)->Mass(); // p mass
fMass[6] = fDB->GetParticle(2212)->Mass()+fDB->GetParticle(2112)->Mass()*2; // p mass
fMass[7] = (fDB->GetParticle(2212)->Mass()+fDB->GetParticle(2112)->Mass()*2)*0.5; // p mass
+ fMass[8] = (fDB->GetParticle(2212)->Mass()*2+fDB->GetParticle(2112)->Mass()*2)*0.5; // p mass
// TOF response
fTOFResponseF = new TF1("fTOFprob","[0]*TMath::Exp(-(x-[1])*(x-[1])/2/[2]/[2])* (x < [1]+[3]*[2]) + (x > [1]+[3]*[2])*[0]*TMath::Exp(-(x-[1]-[3]*[2]*0.5)*[3]/[2])",-7,7);
fPsiRes=999;
}
//________________________________________________________________________
-Float_t AliFlowBayesianPID::GetExpDeDx(const AliESDtrack *t,Int_t iS) const{
+//________________________________________________________________________
+Float_t AliFlowBayesianPID::GetExpDeDx(const AliVTrack *t,Int_t iS) const{
// tuned dE/dx (vs. eta and centrality)
- Double_t ptpc[3];
- t->GetInnerPxPyPz(ptpc);
- Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
+ Float_t momtpc=t->GetTPCmomentum();
Float_t dedxExp=0;
if(iS==0) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kElectron);
else if(iS==5) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kDeuteron);
else if(iS==6) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kTriton);
else if(iS==7) dedxExp = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5;
-
+ else if(iS==8) dedxExp = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[8])*5;
+
Float_t eta = t->Eta();
Float_t etaCorr = 7.98368e-03 - 1.67208e-02 - 1.89776e-01*eta*eta -2.90836e-02*eta*eta + 5.96093e-01*eta*eta*eta*eta + 6.06450e-02*eta*eta*eta*eta - 3.55884e-01*eta*eta*eta*eta*eta*eta;
if(fCurrCentrality < 0){
return dedxExp;
}
//________________________________________________________________________
-Float_t AliFlowBayesianPID::GetExpDeDx(const AliAODTrack *t,Int_t iS) const{
+Float_t AliFlowBayesianPID::GetExpDeDx(const AliVTrack *t,Float_t mass) const{
// tuned dE/dx (vs. eta and centrality)
Float_t momtpc=t->GetTPCmomentum();
- Float_t dedxExp=0;
- if(iS==0) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kElectron);
- else if(iS==1) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kMuon);
- else if(iS==2) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kPion);
- else if(iS==3) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kKaon);
- else if(iS==4) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kProton);
- else if(iS==5) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kDeuteron);
- else if(iS==6) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kTriton);
- else if(iS==7) dedxExp = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5;
+ if(mass < 0.0001) mass = 0.0001;
+
+ Float_t dedxExp = fPIDesd->GetTPCResponse().Bethe(momtpc/mass);
Float_t eta = t->Eta();
Float_t etaCorr = 7.98368e-03 - 1.67208e-02 - 1.89776e-01*eta*eta -2.90836e-02*eta*eta + 5.96093e-01*eta*eta*eta*eta + 6.06450e-02*eta*eta*eta*eta - 3.55884e-01*eta*eta*eta*eta*eta*eta;
+
if(fCurrCentrality < 0){
}
else if(fCurrCentrality < 5) etaCorr += 17E-3;
resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kProton);
dedxExp = GetExpDeDx(t,4);
}
+ else if(iS-1e+9 == 10020){ // d
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kDeuteron);
+ dedxExp = GetExpDeDx(t,5);
+ }
+ else if(iS-1e+9 == 10030){ // t
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kTriton);
+ dedxExp = GetExpDeDx(t,6);
+ }
+ else if(iS-1e+9 == 20030){ // 3He
+ resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5*0.07;
+ dedxExp = GetExpDeDx(t,7);
+ }
+ else if(iS-1e+9 == 20040){ // 4He
+ resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[5])*5*0.07;
+ dedxExp = GetExpDeDx(t,8);
+ }
if(resolutionTPC > -1) dedx = fTPCResponseF->GetRandom()*resolutionTPC + dedxExp;
else dedx = 0;
}
else if(iS==5) resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kDeuteron);
else if(iS==6) resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kTriton);
else if(iS==7) resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5*0.07;
+ else if(iS==8) resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[8])*5*0.07;
if(centr < 0) resolutionTPC *= 0.78;
if(centr < 10) resolutionTPC *= 1.0;
Float_t mismweight = TMath::Max(fgMism->Eval(timeTOF - timeextra),0.0000001) * ((0.5 + 0.05/pt/pt/pt)*(0.75 + 0.23 * (1.3*dx*dx + 0.7*dz*dz))); // mismatch probabilities
fWTofMism = mismfrac*mismweight;
- Double_t inttimes[8];
+ Double_t inttimes[9];
t->GetIntegratedTimes(inttimes);
inttimes[5] = inttimes[0] / p * fMass[5] * TMath::Sqrt(1+p*p/fMass[5]/fMass[5]);
inttimes[6] = inttimes[0] / p * fMass[6] * TMath::Sqrt(1+p*p/fMass[6]/fMass[6]);
inttimes[7] = inttimes[0] / p * fMass[7] * TMath::Sqrt(1+p*p/fMass[7]/fMass[7]);
+ inttimes[8] = inttimes[0] / p * fMass[8] * TMath::Sqrt(1+p*p/fMass[8]/fMass[8]);
for(Int_t iS=0;iS<fgkNspecies;iS++){
Float_t expsigma = fPIDesd->GetTOFResponse().GetExpectedSigma(p, inttimes[iS], fMass[iS]);
resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kPion);
dedxExp = GetExpDeDx(t,2);
}
- else if(iS==321){
- resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kKaon);
- dedxExp = GetExpDeDx(t,3);
- }
- else if(iS==2212){
- resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kProton);
- dedxExp = GetExpDeDx(t,4);
- }
+ else if(iS==321){
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kKaon);
+ dedxExp = GetExpDeDx(t,3);
+ }
+ else if(iS==2212){
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kProton);
+ dedxExp = GetExpDeDx(t,4);
+ }
+ else if(iS-1e+9 == 10020){ // d
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kDeuteron);
+ dedxExp = GetExpDeDx(t,5);
+ }
+ else if(iS-1e+9 == 10030){ // t
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kTriton);
+ dedxExp = GetExpDeDx(t,6);
+ }
+ else if(iS-1e+9 == 20030){ // 3He
+ resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5*0.07;
+ dedxExp = GetExpDeDx(t,7);
+ }
+ else if(iS-1e+9 == 20040){ // 4He
+ resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[5])*5*0.07;
+ dedxExp = GetExpDeDx(t,8);
+ }
if(resolutionTPC > -1)
dedx = fTPCResponseF->GetRandom()*resolutionTPC + dedxExp;
else dedx = 0;
else if(iS==5) resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kDeuteron);
else if(iS==6) resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kTriton);
else if(iS==7) resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5*0.07;
+ else if(iS==8) resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[8])*5*0.07;
if(centr < 0) resolutionTPC *= 0.78;
if(centr < 10) resolutionTPC *= 1.0;
Float_t mismweight = TMath::Max(fgMism->Eval(timeTOF - timeextra),0.0000001) * ((0.5 + 0.05/pt/pt/pt)); // mismatch probabilities
fWTofMism = mismfrac*mismweight;
- Double_t inttimes[8];
+ Double_t inttimes[9];
t->GetIntegratedTimes(inttimes);
inttimes[5] = inttimes[0] / p * fMass[5] * TMath::Sqrt(1+p*p/fMass[5]/fMass[5]);
inttimes[6] = inttimes[0] / p * fMass[6] * TMath::Sqrt(1+p*p/fMass[6]/fMass[6]);
inttimes[7] = inttimes[0] / p * fMass[7] * TMath::Sqrt(1+p*p/fMass[7]/fMass[7]);
+ inttimes[8] = inttimes[0] / p * fMass[8] * TMath::Sqrt(1+p*p/fMass[8]/fMass[8]);
for(Int_t iS=0;iS<fgkNspecies;iS++){
if(iS==0){
fghPriors[5]->SetBinContent(k1,k2,(fC[icentr][TMath::Max(ipt-5,0)][4]*(1-weight) + 0.2*weight)/deutFact);
fghPriors[6]->SetBinContent(k1,k2,(fC[icentr][TMath::Max(ipt-5,0)][4]*(1-weight) + 0.2*weight)/deutFact/deutFact);
fghPriors[7]->SetBinContent(k1,k2,(fC[icentr][TMath::Max(ipt-5,0)][4]*(1-weight) + 0.2*weight)/deutFact/deutFact/deutFact);
- } // end loop on pt
+ fghPriors[8]->SetBinContent(k1,k2,(fC[icentr][TMath::Max(ipt-5,0)][4]*(1-weight) + 0.2*weight)/deutFact/deutFact/deutFact/deutFact);
+ } // end loop on pt
} // end loop on centrality bins
}
Bool_t GetDetANDstatus(Int_t idet) const {if(idet < fgkNdetectors && idet >= 0){return fMaskAND[idet];} else{return kFALSE;} };
Bool_t GetDetORstatus(Int_t idet) const {if(idet < fgkNdetectors && idet >= 0){return fMaskOR[idet];} else{return kFALSE;} };
Bool_t GetCurrentMask(Int_t idet) const {if(idet < fgkNdetectors && idet >= 0){return fMaskCurrent[idet];} else{return kFALSE;} };
- Float_t GetExpDeDx(const AliESDtrack *t,Int_t iS) const;
- Float_t GetExpDeDx(const AliAODTrack *t,Int_t iS) const;
+
+ Float_t GetExpDeDx(const AliVTrack *t,Int_t iS) const;
+ Float_t GetExpDeDx(const AliVTrack *t,Float_t m) const;
// methods for Bayesina Combined PID
void ComputeWeights(const AliESDtrack *t);
void SetPriors();
static const Int_t fgkNdetectors = 2; // Number of detector used for PID
- static const Int_t fgkNspecies = 8;// 0=el, 1=mu, 2=pi, 3=ka, 4=pr, 5=deuteron, 6=triton, 7=He3
+ static const Int_t fgkNspecies = 9;// 0=el, 1=mu, 2=pi, 3=ka, 4=pr, 5=deuteron, 6=triton, 7=He3
static TH2D* fghPriors[fgkNspecies]; // histo with priors (hardcoded)
static TSpline3 *fgMism; // function for mismatch
static TH1D *fgHtofChannelDist; // channel distance from IP
- ClassDef(AliFlowBayesianPID, 7); // example of analysis
+ ClassDef(AliFlowBayesianPID, 8); // example of analysis
};
#endif
protected:
Int_t fNDaughters; // number of daughters (5 max)
Int_t fDaughter[5]; // fID of daughter, points back to ESD track
- AliFlowTrack *fTrack[5]; // pointer to daughter in FlowEvent
-
+ AliFlowTrack *fTrack[5]; //! pointer to daughter in FlowEvent
ClassDef(AliFlowCandidateTrack, 2);
};
void AddSpecies(const char *name,Int_t nXbin,const Double_t *xbin,Int_t nYbin,const Double_t *ybin);
- const char *GetSpeciesName(Int_t species){return GetQA(species)->GetName();};
+ const char *GetSpeciesName(Int_t species){if(species<GetNspecies()) return GetQA(species*GetNhistos()/GetNspecies())->GetName();else return "";};
Int_t Add(const AliFlowVZEROQA *oth);
char nameHisto[200];
char title[300];
+ char title2[300];
for(Int_t i=0; i < ncomb;i++){
snprintf(nameHisto,200,"%s_%s_%i",GetName(),name,i);
snprintf(title,300,"%s",name);
Int_t ncombTemp = i;
for(Int_t j=0;j < GetNvar();j++){
Int_t ibin = ncombTemp%(*fNbinVar)[j];
- snprintf(title,300,"%s_%04.1f<%s<%04.1f",title,(*fXmin)[j] + ((*fXmax)[j]-(*fXmin)[j])/(*fNbinVar)[j]*ibin,fNameVar->At(j)->GetName(),(*fXmin)[j] + ((*fXmax)[j]-(*fXmin)[j])/(*fNbinVar)[j]*(ibin+1));
+ snprintf(title2,300,"%s",title);
+ snprintf(title,300,"%s_%04.1f<%s<%04.1f",title2,(*fXmin)[j] + ((*fXmax)[j]-(*fXmin)[j])/(*fNbinVar)[j]*ibin,fNameVar->At(j)->GetName(),(*fXmin)[j] + ((*fXmax)[j]-(*fXmin)[j])/(*fNbinVar)[j]*(ibin+1));
ncombTemp /= (*fNbinVar)[j];
}
return GetV2(species);
}
+TProfile *AliFlowVZEROResults::GetV2reweight(Int_t species,Float_t xMin[],Float_t xMax[],Int_t varRW,Float_t rw[]) const{
+ if(GetNvar()){
+ char title[300];
+ char title2[300];
+ Int_t ncomb = 1;
+ for(Int_t i=0;i < GetNvar();i++){
+ ncomb *= (*fNbinVar)[i];
+ }
+
+ TProfile *htemplate = GetV2(species*ncomb);
+ TProfile *temp = new TProfile(*htemplate);
+ temp->SetName("histo");
+ temp->Reset();
+ snprintf(title,300,"%i",species);
+ for(Int_t i=0;i < GetNvar();i++){
+ Int_t imin = GetBin(i,xMin[i]);
+ if(imin < 0) imin = 0;
+ else if(imin >= (*fNbinVar)[i]) imin = (*fNbinVar)[i]-1;
+ Int_t imax = GetBin(i,xMax[i]);
+ if(imax < imin) imax = imin;
+ else if(imax >= (*fNbinVar)[i]) imax = (*fNbinVar)[i]-1;
+ snprintf(title2,300,"%s",title);
+ snprintf(title,300,"%s_%04.1f<%s<%04.1f",title2,
+ (*fXmin)[i] + ((*fXmax)[i]-(*fXmin)[i])/(*fNbinVar)[i]*imin,
+ fNameVar->At(i)->GetName(),
+ (*fXmin)[i] + ((*fXmax)[i]-(*fXmin)[i])/(*fNbinVar)[i]*(imax+1));
+ }
+ temp->SetTitle(title);
+
+ for(Int_t i=species*ncomb;i < (species+1)*ncomb;i++){
+ Bool_t kGood = kTRUE;
+
+ Int_t currentBin=0;
+
+ Int_t ncombTemp = i;
+ for(Int_t j=0;j < GetNvar();j++){
+ Int_t imin = GetBin(j,xMin[j]);
+ if(imin < 0) imin = 0;
+ else if(imin >= (*fNbinVar)[j]) imin = (*fNbinVar)[j]-1;
+ Int_t imax = GetBin(j,xMax[j]);
+ if(imax < imin) imax = imin;
+ else if(imax >= (*fNbinVar)[j]) imax = (*fNbinVar)[j]-1;
+
+ Int_t ibin = ncombTemp%(*fNbinVar)[j];
+ ncombTemp /= (*fNbinVar)[j];
+
+ if(j == varRW) currentBin = ibin;
+
+ if(ibin < imin || ibin > imax){
+ kGood = kFALSE;
+ j = GetNvar();
+ }
+ }
+
+ if(kGood) temp->Add(GetV2(i),rw[currentBin]);
+ }
+ return temp;
+ }
+
+ return GetV2(species);
+}
+
Long64_t AliFlowVZEROResults::Merge(TCollection* list){
Long64_t res=0;
TProfile *GetV2(Int_t histo) const {return ((TProfile *) fV2->At(histo));};
TProfile *GetV2(Int_t species,Float_t x[]) const;
TProfile *GetV2(Int_t species,Float_t xMin[],Float_t xMax[]) const;
+ TProfile *GetV2reweight(Int_t species,Float_t xMin[],Float_t xMax[],Int_t varRW,Float_t rw[]) const;
void DirectFill(Int_t histo,Float_t pt,Float_t v2){GetV2(histo)->Fill(pt,v2);};
void Fill(Int_t species,Float_t pt,Float_t v2,Float_t x[]);
void AddSpecies(const char *name,Int_t nXbin,const Double_t *bin);
- const char *GetSpeciesName(Int_t species){return GetV2(species)->GetName();};
+ const char *GetSpeciesName(Int_t species){if(species < GetNspecies()) return GetV2(species*GetNhistos()/GetNspecies())->GetName();else return "";};
Int_t Add(const AliFlowVZEROResults *oth);
-AliAnalysisTask *AddTaskVZERO(Bool_t ismc=kFALSE,Bool_t kV2=kTRUE,Bool_t kV3=kTRUE,Bool_t qa=kTRUE){
+AliAnalysisTask *AddTaskVZERO(Bool_t ismc=kFALSE,Bool_t kV2=kTRUE,Bool_t kV3=kTRUE,Bool_t qa=kTRUE,Bool_t modulationdEdx=kFALSE,Bool_t globalTrack=kFALSE){
//get the current analysis manager
AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
if(ismc) task->SetMC();
if(qa) task->SetQA();
+ task->SetFillDCAinfo(globalTrack); // 0 = TPC only track, 1 = global tracks
+ task->SetModulationDEDx(modulationdEdx);
+
mgr->AddTask(task);
//Attach input to my tasks
-extractFlowVZERO(Int_t icentr,Int_t spec,Int_t arm=2,Bool_t isMC=kFALSE){
+TFile *fo = NULL;
+TF1 *fPsi = NULL;
+
+Bool_t kLib = kFALSE;
+Bool_t kVZEROrescorr = kTRUE; // VZERO resolution corrections
+Bool_t kNUAcorr = kTRUE; // NUA corrections
+Bool_t kNUOcorr = kFALSE; // Not Uniform Occupancy (TOF) corrections
+Bool_t kPIDcorr = kFALSE; // PID corrections
+Bool_t kCleanMemory = kTRUE; // if kFALSE take care of the large numbers of TCanvases
+
+TH1D *hNUO[8]; // NUO corrections
+
+void LoadLib();
+void extractFlowVZERO(Int_t icentr,Int_t arm=2,Float_t pTh=0.8,Bool_t isMC=kFALSE,Int_t addbin=0);
+TProfile *extractFlowVZEROsingle(Int_t icentr,Int_t spec,Int_t arm=2,Bool_t isMC=kFALSE,Float_t pTh=0.9,Int_t addbin=0,const char *nameSp="",Float_t detMin=0,Float_t detMax=1,Int_t chMin=-1,Int_t chMax=1);
+void compareTPCTOF(Int_t icentr,Int_t spec,Int_t arm=2,Float_t pTh=0.9,Int_t addbin=0);
+void plotQApid(Int_t ic,Float_t pt,Int_t addbin=0);
+
+void LoadLib(){
+ if(! kLib){
+ gSystem->Load("libVMC.so");
+ gSystem->Load("libPhysics.so");
+ gSystem->Load("libTree.so");
+ gSystem->Load("libSTEERBase.so");
+ gSystem->Load("libANALYSIS.so");
+ gSystem->Load("libAOD.so");
+ gSystem->Load("libESD.so");
+ gSystem->Load("libANALYSIS.so");
+ gSystem->Load("libANALYSISalice.so");
+ gSystem->Load("libCORRFW.so");
+ gSystem->Load("libNetx.so");
+ gSystem->Load("libPWGflowBase.so");
+ gSystem->Load("libPWGflowTasks.so");
+
+ kLib = kTRUE;
+ }
+}
+
+void extractFlowVZERO(Int_t icentr,Int_t arm,Float_t pTh,Bool_t isMC,Int_t addbin){
+ LoadLib();
+
+ new TCanvas();
+
+ Int_t cMin[] = {0,5,10,20,30,40,50,60,70};
+ Int_t cMax[] = {5,10,20,30,40,50,60,70,80};
+
+ if(kNUOcorr){ // Compute correction for NUO in TOF
+ compareTPCTOF(icentr,0,arm,pTh,addbin);
+// compareTPCTOF(icentr,1,arm,pTh,addbin);
+// compareTPCTOF(icentr,2,arm,pTh,addbin);
+// compareTPCTOF(icentr,3,arm,pTh,addbin);
+// compareTPCTOF(icentr,4,arm,pTh,addbin);
+ }
+
+ TProfile *pAll;
+ pAll=extractFlowVZEROsingle(icentr,0,arm,0,pTh,addbin,"all",0,1);
+ pAll->SetMarkerStyle(24);
+ TProfile *pPiTOF,*pPiTPC,*pPiTPC2;
+ pPiTOF=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTOF",1,1);
+ pPiTPC=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTPC",0,0);
+ pPiTPC2=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTPC2",2,2);
+ pPiTPC->Add(pPiTPC2);
+
+ TH1D *hPi = pPiTOF->ProjectionX("hPi");
+ hPi->SetLineColor(4);
+ hPi->SetMarkerColor(4);
+ hPi->SetMarkerStyle(20);
+ for(Int_t i=1;i <=hPi->GetNbinsX();i++){
+ Float_t x = hPi->GetBinCenter(i);
+ if(x < 0.2){
+ hPi->SetBinContent(i,0);
+ hPi->SetBinError(i,0);
+ }
+ else if(x < 0.5){
+ hPi->SetBinContent(i,pPiTPC->GetBinContent(i));
+ hPi->SetBinError(i,pPiTPC->GetBinError(i));
+ }
+ else{
+ if(kNUOcorr){
+ hPi->SetBinContent(i,pPiTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
+ hPi->SetBinError(i,pPiTOF->GetBinError(i));
+ }
+ else{
+ hPi->SetBinContent(i,pPiTOF->GetBinContent(i));
+ hPi->SetBinError(i,pPiTOF->GetBinError(i));
+ }
+ }
+ }
+
+ TProfile *pElTOF,*pElTPC,*pElTPC2;
+ pElTOF=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTOF",1,1);
+ pElTPC=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTPC",0,0);
+ pElTPC2=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTPC2",2,2);
+ pElTPC->Add(pElTPC2);
+
+ TH1D *hEl = pElTOF->ProjectionX("hEl");
+ hEl->SetLineColor(6);
+ hEl->SetMarkerColor(6);
+ hEl->SetMarkerStyle(20);
+ for(Int_t i=1;i <=hEl->GetNbinsX();i++){
+ Float_t x = hEl->GetBinCenter(i);
+ if(x < 0.2){
+ hEl->SetBinContent(i,0);
+ hEl->SetBinError(i,0);
+ }
+ else if(x < 0.3){
+ hEl->SetBinContent(i,pElTPC->GetBinContent(i));
+ hEl->SetBinError(i,pElTPC->GetBinError(i));
+ }
+ else{
+ if(kNUOcorr){
+ hEl->SetBinContent(i,pElTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
+ hEl->SetBinError(i,pElTOF->GetBinError(i));
+ }
+ else{
+ hEl->SetBinContent(i,pElTOF->GetBinContent(i));
+ hEl->SetBinError(i,pElTOF->GetBinError(i));
+ }
+ }
+ }
+
+ TProfile *pKTOF,*pKTPC,*pKTPC2;
+ pKTOF=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTOF",1,1);
+ pKTPC=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTPC",0,0);
+ pKTPC2=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTPC2",2,2);
+ pKTPC->Add(pKTPC2);
+
+ TH1D *hK = pKTOF->ProjectionX("hKa");
+ hK->SetLineColor(1);
+ hK->SetMarkerColor(1);
+ hK->SetMarkerStyle(21);
+ for(Int_t i=1;i <=hK->GetNbinsX();i++){
+ Float_t x = hK->GetBinCenter(i);
+ if(x < 0.25){
+ hK->SetBinContent(i,0);
+ hK->SetBinError(i,0);
+ }
+ else if(x < 0.45){
+ hK->SetBinContent(i,pKTPC->GetBinContent(i));
+ hK->SetBinError(i,pKTPC->GetBinError(i));
+ }
+ else{
+ if(kNUOcorr){
+ hK->SetBinContent(i,pKTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
+ hK->SetBinError(i,pKTOF->GetBinError(i));
+ }
+ else{
+ hK->SetBinContent(i,pKTOF->GetBinContent(i));
+ hK->SetBinError(i,pKTOF->GetBinError(i));
+ }
+ }
+ }
+
+ TProfile *pPrTOF,*pPrTOF2,*pPrTPC,*pPrTPC2;
+ pPrTOF=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTOF",1,1,-1,-1);
+ pPrTOF2=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTOF2",1,1,-1,1);
+ pPrTPC=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTPC",0,0,-1,-1);
+ pPrTPC2=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTPC2",2,2,-1,-1);
+ pPrTPC->Add(pPrTPC2);
+
+ TH1D *hPr = pPrTOF->ProjectionX("hPr");
+ hPr->SetLineColor(2);
+ hPr->SetMarkerColor(2);
+ hPr->SetMarkerStyle(22);
+ for(Int_t i=1;i <=hPr->GetNbinsX();i++){
+ Float_t x = hPr->GetBinCenter(i);
+ if(x < 0.3){
+ hPr->SetBinContent(i,0);
+ hPr->SetBinError(i,0);
+ }
+ else if(x < 1.0){
+ hPr->SetBinContent(i,pPrTPC->GetBinContent(i));
+ hPr->SetBinError(i,pPrTPC->GetBinError(i));
+ }
+ else{
+ if(x < 3){
+ if(kNUOcorr){
+ hPr->SetBinContent(i,pPrTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
+ hPr->SetBinError(i,pPrTOF->GetBinError(i));
+ }
+ else{
+ hPr->SetBinContent(i,pPrTOF->GetBinContent(i));
+ hPr->SetBinError(i,pPrTOF->GetBinError(i));
+ }
+ }
+ else{
+ if(kNUOcorr){
+ hPr->SetBinContent(i,pPrTOF2->GetBinContent(i) + hNUO[0]->GetBinContent(i));
+ hPr->SetBinError(i,pPrTOF2->GetBinError(i));
+ }
+ else{
+ hPr->SetBinContent(i,pPrTOF2->GetBinContent(i));
+ hPr->SetBinError(i,pPrTOF2->GetBinError(i));
+ }
+ }
+ }
+ }
+
+ pAll->Draw();
+ hPi->Draw("SAME");
+ hK->Draw("SAME");
+ hPr->Draw("SAME");
+
+
+ char name[100];
+
+ // PID correction
+ if(kPIDcorr){
+ TFile *fPidTOF = new TFile("$ALICE_ROOT/PWGCF/FLOW/other/BayesianPIDcontTPCTOF.root");
+ TFile *fPidTPC = new TFile("$ALICE_ROOT/PWGCF/FLOW/other/BayesianPIDcontTPC.root");
+ // pi histos
+ sprintf(name,"Pi_IDas_Picentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPiPi=(TH1D *) fPidTOF->Get(name);
+ sprintf(name,"Pi_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPiEl=(TH1D *) fPidTOF->Get(name);
+ sprintf(name,"Pi_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPiKa=(TH1D *) fPidTOF->Get(name);
+ sprintf(name,"Pi_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPiPr=(TH1D *) fPidTOF->Get(name);
+ TH1D *hPidAll = new TH1D(*hPidPiPi);
+ hPidAll->Add(hPidPiKa);
+ hPidAll->Add(hPidPiPr);
+ hPidAll->Add(hPidPiEl);
+ hPidPiPi->Divide(hPidAll);
+ hPidPiKa->Divide(hPidAll);
+ hPidPiPr->Divide(hPidAll);
+ hPidPiEl->Divide(hPidAll);
+
+ sprintf(name,"Pi_IDas_Picentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPiPiTPC=(TH1D *) fPidTPC->Get(name);
+ sprintf(name,"Pi_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPiElTPC=(TH1D *) fPidTPC->Get(name);
+ sprintf(name,"Pi_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPiKaTPC=(TH1D *) fPidTPC->Get(name);
+ sprintf(name,"Pi_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPiPrTPC=(TH1D *) fPidTPC->Get(name);
+ hPidAll->Reset();
+ hPidAll->Add(hPidPiPiTPC);
+ hPidAll->Add(hPidPiKaTPC);
+ hPidAll->Add(hPidPiPrTPC);
+ hPidAll->Add(hPidPiElTPC);
+ hPidPiPiTPC->Divide(hPidAll);
+ hPidPiKaTPC->Divide(hPidAll);
+ hPidPiPrTPC->Divide(hPidAll);
+ hPidPiElTPC->Divide(hPidAll);
+
+ // K histos
+ sprintf(name,"Ka_IDas_Picentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidKaPi=(TH1D *) fPidTOF->Get(name);
+ sprintf(name,"Ka_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidKaEl=(TH1D *) fPidTOF->Get(name);
+ sprintf(name,"Ka_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidKaKa=(TH1D *) fPidTOF->Get(name);
+ sprintf(name,"Ka_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidKaPr=(TH1D *) fPidTOF->Get(name);
+ hPidAll->Reset();
+ hPidAll->Add(hPidKaPi);
+ hPidAll->Add(hPidKaKa);
+ hPidAll->Add(hPidKaPr);
+ hPidAll->Add(hPidKaEl);
+ hPidKaPi->Divide(hPidAll);
+ hPidKaKa->Divide(hPidAll);
+ hPidKaPr->Divide(hPidAll);
+ hPidKaEl->Divide(hPidAll);
+
+ sprintf(name,"Ka_IDas_Picentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidKaPiTPC=(TH1D *) fPidTPC->Get(name);
+ sprintf(name,"Ka_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidKaElTPC=(TH1D *) fPidTPC->Get(name);
+ sprintf(name,"Ka_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidKaKaTPC=(TH1D *) fPidTPC->Get(name);
+ sprintf(name,"Ka_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidKaPrTPC=(TH1D *) fPidTPC->Get(name);
+ hPidAll->Reset();
+ hPidAll->Add(hPidKaPiTPC);
+ hPidAll->Add(hPidKaKaTPC);
+ hPidAll->Add(hPidKaPrTPC);
+ hPidAll->Add(hPidKaElTPC);
+ hPidKaPiTPC->Divide(hPidAll);
+ hPidKaKaTPC->Divide(hPidAll);
+ hPidKaPrTPC->Divide(hPidAll);
+ hPidKaElTPC->Divide(hPidAll);
+
+ // pr histos
+ sprintf(name,"Pr_IDas_Picentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPrPi=(TH1D *) fPidTOF->Get(name);
+ sprintf(name,"Pr_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPrEl=(TH1D *) fPidTOF->Get(name);
+ sprintf(name,"Pr_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPrKa=(TH1D *) fPidTOF->Get(name);
+ sprintf(name,"Pr_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPrPr=(TH1D *) fPidTOF->Get(name);
+ hPidAll->Reset();
+ hPidAll->Add(hPidPrPi);
+ hPidAll->Add(hPidPrKa);
+ hPidAll->Add(hPidPrPr);
+ hPidAll->Add(hPidPrEl);
+ hPidPrPi->Divide(hPidAll);
+ hPidPrKa->Divide(hPidAll);
+ hPidPrPr->Divide(hPidAll);
+ hPidPrEl->Divide(hPidAll);
+
+ sprintf(name,"Pr_IDas_Picentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPrPiTPC=(TH1D *) fPidTPC->Get(name);
+ sprintf(name,"Pr_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPrElTPC=(TH1D *) fPidTPC->Get(name);
+ sprintf(name,"Pr_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPrKaTPC=(TH1D *) fPidTPC->Get(name);
+ sprintf(name,"Pr_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
+ TH1D *hPidPrPrTPC=(TH1D *) fPidTPC->Get(name);
+ hPidAll->Reset();
+ hPidAll->Add(hPidPrPiTPC);
+ hPidAll->Add(hPidPrKaTPC);
+ hPidAll->Add(hPidPrPrTPC);
+ hPidAll->Add(hPidPrElTPC);
+ hPidPrPiTPC->Divide(hPidAll);
+ hPidPrKaTPC->Divide(hPidAll);
+ hPidPrPrTPC->Divide(hPidAll);
+ hPidPrElTPC->Divide(hPidAll);
+
+ for(Int_t k=1;k <=hPi->GetNbinsX();k++){
+ Float_t pt = hPi->GetBinCenter(k);
+ Float_t xPi=hPi->GetBinContent(k)*hPidPiPi->Interpolate(pt) + hK->GetBinContent(k)*hPidPiKa->Interpolate(pt) + hPr->GetBinContent(k)*hPidPiPr->Interpolate(pt) + hEl->GetBinContent(k)*hPidPiEl->Interpolate(pt);
+ if(pt < 0.5)
+ xPi=hPi->GetBinContent(k)*hPidPiPiTPC->Interpolate(pt) + hK->GetBinContent(k)*hPidPiKaTPC->Interpolate(pt) + hPr->GetBinContent(k)*hPidPiPrTPC->Interpolate(pt) + hEl->GetBinContent(k)*hPidPiElTPC->Interpolate(pt);
+ Float_t xKa=hPi->GetBinContent(k)*hPidKaPi->Interpolate(pt) + hK->GetBinContent(k)*hPidKaKa->Interpolate(pt) + hPr->GetBinContent(k)*hPidKaPr->Interpolate(pt) + hEl->GetBinContent(k)*hPidKaEl->Interpolate(pt);
+ if(pt < 0.45)
+ xKa=hPi->GetBinContent(k)*hPidKaPiTPC->Interpolate(pt) + hK->GetBinContent(k)*hPidKaKaTPC->Interpolate(pt) + hPr->GetBinContent(k)*hPidKaPrTPC->Interpolate(pt) + hEl->GetBinContent(k)*hPidKaElTPC->Interpolate(pt);
+ Float_t xPr=hPi->GetBinContent(k)*hPidPrPi->Interpolate(pt) + hK->GetBinContent(k)*hPidPrKa->Interpolate(pt) + hPr->GetBinContent(k)*hPidPrPr->Interpolate(pt) + hEl->GetBinContent(k)*hPidPrEl->Interpolate(pt);
+ if(pt < 1)
+ xPr=hPi->GetBinContent(k)*hPidPrPiTPC->Interpolate(pt) + hK->GetBinContent(k)*hPidPrKaTPC->Interpolate(pt) + hPr->GetBinContent(k)*hPidPrPrTPC->Interpolate(pt) + hEl->GetBinContent(k)*hPidPrElTPC->Interpolate(pt);
+
+ hPi->SetBinContent(k,hPi->GetBinContent(k)*2 - xPi);
+ hK->SetBinContent(k,hK->GetBinContent(k)*2 - xKa);
+ hPr->SetBinContent(k,hPr->GetBinContent(k)*2 - xPr);
+ }
+ }
+
+ // write output
+ snprintf(name,100,"results%03i-%03iv%i_pTh%3.1f.root",cMin[icentr],cMax[icentr+addbin],arm,pTh);
+ TFile *fout = new TFile(name,"RECREATE");
+ pAll->ProjectionX()->Write();
+ hPi->Write();
+ hK->Write();
+ hPr->Write();
+ if(isMC){
+ TH1D *pTmp = extractFlowVZEROsingle(icentr,0,arm,1,pTh,addbin,"allMC",1,1,-1,1)->ProjectionX();
+ pTmp->SetLineColor(6);
+ pTmp->SetMarkerColor(6);
+ pTmp->SetMarkerStyle(24);
+ pTmp->Write();
+ pTmp = extractFlowVZEROsingle(icentr,1,arm,1,pTh,addbin,"piMC",1,1,-1,1)->ProjectionX();
+ pTmp->SetLineColor(4);
+ pTmp->SetMarkerColor(4);
+ pTmp->SetMarkerStyle(24);
+ pTmp->Write();
+ pTmp = extractFlowVZEROsingle(icentr,2,arm,1,pTh,addbin,"kaMC",1,1,-1,1)->ProjectionX();
+ pTmp->SetLineColor(1);
+ pTmp->SetMarkerColor(1);
+ pTmp->SetMarkerStyle(25);
+ pTmp->Write();
+ pTmp = extractFlowVZEROsingle(icentr,3,arm,1,pTh,addbin,"prMC",1,1,-1,-1)->ProjectionX();
+ pTmp->SetLineColor(2);
+ pTmp->SetMarkerColor(2);
+ pTmp->SetMarkerStyle(26);
+ pTmp->Write();
+ }
+ fout->Close();
+}
+
+TProfile *extractFlowVZEROsingle(Int_t icentr,Int_t spec,Int_t arm,Bool_t isMC,Float_t pTh,Int_t addbin,const char *nameSp,Float_t detMin,Float_t detMax,Int_t chMin,Int_t chMax){
+ LoadLib();
+
+ pTh += 0.00001;
+
// NUA correction currently are missing
char name[100];
- snprintf(name,100,"outVZEROv%i.root",arm);
- TFile *fo = new TFile(name);
+ char stringa[200];
+
+ snprintf(name,100,"AnalysisResults.root");
+ if(!fo) fo = new TFile(name);
snprintf(name,100,"contVZEROv%i",arm);
TList *cont = (TList *) fo->Get(name);
- Float_t xMin[5] = {icentr,-1,0,-10,0};
- Float_t xMax[5] = {icentr,1,1,10,1.5};
+ cont->ls();
+
+ Float_t xMin[5] = {icentr/*centrality bin*/,chMin/*charge*/,pTh/*prob*/,-TMath::Pi()/arm/*Psi*/,detMin/*PID mask*/};
+ Float_t xMax[5] = {icentr+addbin,chMax,1,TMath::Pi()/arm,detMax};
cont->ls();
TProfile *p2 = cont->At(3);
TProfile *p3 = cont->At(4);
+ TH2F *hPsi2DA = cont->At(5);
+ TH2F *hPsi2DC = cont->At(6);
+
+ TH1D *hPsiA = hPsi2DA->ProjectionY("PsiA",icentr+1,icentr+addbin+1);
+ TH1D *hPsiC = hPsi2DC->ProjectionY("PsiC",icentr+1,icentr+addbin+1);
+
+ if(!fPsi) fPsi = new TF1("fPsi","pol0",-TMath::Pi()/arm,TMath::Pi()/arm);
+ hPsiA->Fit(fPsi,"0");
+ Float_t offsetA = fPsi->GetParameter(0);
+ hPsiC->Fit(fPsi,"0");
+ Float_t offsetC = fPsi->GetParameter(0);
+
+ Int_t nbinPsi = hPsiA->GetNbinsX();
+
+ Float_t *NUAcorrA = new Float_t[nbinPsi];
+ Float_t *NUAcorrC = new Float_t[nbinPsi];
+
+ for(Int_t i=0;i < nbinPsi;i++){
+ NUAcorrA[i] = offsetA/(hPsiA->GetBinContent(i+1));
+ NUAcorrC[i] = offsetC/(hPsiC->GetBinContent(i+1));
+ }
+
Float_t res1=0,res2=0,res3=0;
Float_t eres1=0,eres2=0,eres3=0;
- Int_t i = icentr;
- 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);
+ for(Int_t i = icentr; i <= icentr+addbin;i++){
+ 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(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);
}
-
+
res1 /= eres1;
res2 /= eres2;
res3 /= eres3;
+ eres1 = sqrt(1./eres1);
+ eres2 = sqrt(1./eres2);
+ eres3 = sqrt(1./eres3);
+
AliFlowVZEROResults *a = (AliFlowVZEROResults *) cont->At(0);
AliFlowVZEROResults *b = (AliFlowVZEROResults *) cont->At(1);
- TProfile *pp = a->GetV2(spec,xMin,xMax);
- TProfile *pp2 = b->GetV2(spec,xMin,xMax);
-
+ TProfile *pp,*pp2;
+ if(kNUAcorr){ // with NUA corrections
+ pp = a->GetV2reweight(spec,xMin,xMax,3,NUAcorrA);
+ pp2 = b->GetV2reweight(spec,xMin,xMax,3,NUAcorrC);
+ }
+ else{
+ pp = a->GetV2(spec,xMin,xMax);
+ pp2 = b->GetV2(spec,xMin,xMax);
+ }
Float_t scaling = sqrt(res1*res3/res2);
- pp->Scale(1./scaling);
-
- printf("resolution V0A = %f\n",scaling);
+ if(kVZEROrescorr){
+ pp->Scale(1./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;
+ printf("resolution V0A = %f +/- %f\n",scaling,err1_2);
scaling = sqrt(res2*res3/res1);
err2_2 /= scaling*scaling;
- pp2->Scale(1./scaling);
- printf("resolution V0C =%f\n",scaling);
+ if(kVZEROrescorr){
+ pp2->Scale(1./scaling);
+ }
+ printf("resolution V0C =%f +/- %f\n",scaling,err2_2);
pp->SetName("V0A");
pp2->SetName("V0C");
- pp->Draw();
- pp2->Draw("SAME");
+ if(! kCleanMemory){
+ new TCanvas();
+ pp->Draw();
+ pp2->Draw("SAME");
+ }
+
+ TProfile *pData = new TProfile(*pp);
+ pData->Add(pp2);
+ snprintf(stringa,100,"%sData",nameSp);
+ pData->SetName(stringa);
- if(arm == 2 && isMC){
- snprintf(name,100,"outVZEROmc.root");
- fo = new TFile(name);
+ TProfile *pMc = NULL;
+
+ TProfile *ppMC;
+ TProfile *ppMC2;
+
+ if(isMC){
snprintf(name,100,"contVZEROmc");
cont = (TList *) fo->Get(name);
- AliFlowVZEROResults *c = (AliFlowVZEROResults *) cont->At(0);
- c->GetV2(spec,xMin,xMax)->Draw("SAME");
+ cont->ls();
+ if(arm == 2){
+ AliFlowVZEROResults *c = (AliFlowVZEROResults *) cont->At(0);
+ if(! kCleanMemory) c->GetV2(spec,xMin,xMax)->Draw("SAME");
+ }
+ AliFlowVZEROResults *cA;
+ if(fo->Get("contVZEROv2")) cA = (AliFlowVZEROResults *) cont->At(1+2*(arm==3));
+ else cA = (AliFlowVZEROResults *) cont->At(0);
+ AliFlowVZEROResults *cC;
+ if(fo->Get("contVZEROv2")) cC = (AliFlowVZEROResults *) cont->At(2+2*(arm==3));
+ else cC = (AliFlowVZEROResults *) cont->At(1);
+
+ TProfile *p1mc = cont->At(5+3*(arm==3));
+ TProfile *p2mc = cont->At(6+3*(arm==3));
+ TProfile *p3mc = cont->At(7+3*(arm==3));
+
+ Float_t resMC1=0,resMC2=0,resMC3=0;
+ Float_t eresMC1=0,eresMC2=0,eresMC3=0;
+
+ for(Int_t i = icentr; i <= icentr+addbin;i++){
+ if(p1mc->GetBinError(i+1)){
+ eresMC1 += 1./p1mc->GetBinError(i+1)/p1mc->GetBinError(i+1);
+ resMC1 += p1mc->GetBinContent(i+1)/p1mc->GetBinError(i+1)/p1mc->GetBinError(i+1);
+ }
+ if(p2mc->GetBinError(i+1)){
+ eresMC2 += 1./p2mc->GetBinError(i+1)/p2mc->GetBinError(i+1);
+ resMC2 += p2mc->GetBinContent(i+1)/p2mc->GetBinError(i+1)/p2mc->GetBinError(i+1);
+ }
+ if(p3mc->GetBinError(i+1)){
+ eresMC3 += 1./p3mc->GetBinError(i+1)/p3mc->GetBinError(i+1);
+ resMC3 += p3mc->GetBinContent(i+1)/p3mc->GetBinError(i+1)/p3mc->GetBinError(i+1);
+ }
+ }
+
+ resMC1 /= eresMC1;
+ resMC2 /= eresMC2;
+ resMC3 /= eresMC3;
+
+ eresMC1 = sqrt(1./eresMC1);
+ eresMC2 = sqrt(1./eresMC2);
+ eresMC3 = sqrt(1./eresMC3);
+
+ ppMC = cA->GetV2(spec,xMin,xMax);
+ ppMC2 = cC->GetV2(spec,xMin,xMax);
+
+ scaling = sqrt(resMC1*resMC3/resMC2);
+ ppMC->Scale(1./scaling);
+
+ err1_2 = eresMC1*eresMC1/resMC1/resMC1/4 +
+ eresMC2*eresMC2/resMC2/resMC2/4 +
+ eresMC3*eresMC3/resMC3/resMC3/4;
+ err2_2 = err1_2;
+ err1_2 /= scaling*scaling;
+ printf("resolution V0A (MC) = %f +/- %f\n",scaling,err1_2);
+ scaling = sqrt(resMC2*resMC3/resMC1);
+ err2_2 /= scaling*scaling;
+ ppMC2->Scale(1./scaling);
+ printf("resolution V0C (MC) =%f +/- %f\n",scaling,err2_2);
+
+ ppMC->SetName("V0Amc");
+ ppMC2->SetName("V0Cmc");
+
+ if(! kCleanMemory){
+ ppMC->Draw("SAME");
+ ppMC2->Draw("SAME");
+ }
+
+ pMc = new TProfile(*ppMC);
+ pMc->Add(ppMC2);
+ snprintf(stringa,100,"%sMC",nameSp);
+ pMc->SetName(stringa);
+ pMc->SetLineColor(2);
}
+
+ if(! kCleanMemory){
+ new TCanvas();
+ pData->Draw();
+ }
+ if(pMc && !kCleanMemory){
+ pMc->Draw("SAME");
+ TH1D *hData = pData->ProjectionX();
+ TH1D *hMc = pMc->ProjectionX();
+ hData->Divide(hMc);
+ new TCanvas();
+ hData->Draw();
+ }
+
+ delete[] NUAcorrA;
+ delete[] NUAcorrC;
+
+ if(kCleanMemory){
+ if(pp) delete pp;
+ if(pp2) delete pp2;
+ if(ppMC) delete ppMC;
+ if(ppMC2) delete ppMC2;
+ }
+
+ delete cont;
+
+ if(isMC) return pMc;
+ return pData;
+}
+void compareTPCTOF(Int_t icentr,Int_t spec,Int_t arm,Float_t pTh,Int_t addbin){
+ LoadLib();
+ Float_t ptMaxFit[8] = {20,0.7,0.55,1.2,2,2,2,2};
+
+ TProfile *pTPC = extractFlowVZEROsingle(icentr,spec,arm,0,pTh,addbin,"TPC",0,0,-1,-1+2*(spec!=3)); // TPC && !TOF (TPC PID)
+ TProfile *pTPC2 = extractFlowVZEROsingle(icentr,spec,arm,0,pTh,addbin,"TPCextra",2,2,-1,-1+2*(spec!=3)); //TPC && TOF (TPC PID)
+
+ TProfile *pTOF = extractFlowVZEROsingle(icentr,spec,arm,0,pTh,addbin,"TOF",1,1,-1,-1+2*(spec!=3)); //TPC && TOF (TPC&TOF PID)
+
+ if(spec > 0) pTPC->Add(pTPC2);
+ else pTPC->Add(pTOF);
+
+ char name[100];
+ snprintf(name,100,"NUO_%i",spec);
+
+ hNUO[spec] = pTPC->ProjectionX(name);
+ hNUO[spec]->Add(pTOF,-1);
+
+ if(spec && hNUO[0]) hNUO[spec]->Add(hNUO[0],-1);
+
+ if(! kCleanMemory){
+ new TCanvas();
+ pTPC->Draw();
+ pTOF->Draw("SAME");
+
+ new TCanvas();
+ pTPC->SetLineColor(1);
+ pTOF->SetLineColor(2);
+ }
+
+ hNUO[spec]->Draw();
+
+ if(kCleanMemory){
+ if(pTPC) delete pTPC;
+ if(pTPC2) delete pTPC2;
+ if(pTOF) delete pTOF;
+ }
+}
+
+void plotQApid(Int_t ic,Float_t pt,Int_t addbin){
+ char name[100];
+ char stringa[200];
+ LoadLib();
+
+ snprintf(name,100,"AnalysisResults.root");
+ if(!fo) fo = new TFile(name);
+ snprintf(name,100,"contVZEROv2");
+ TList *cont = (TList *) fo->Get(name);
+ AliFlowVZEROResults *pidqa = cont->FindObject("qaPID");
+ Float_t xval[2] = {2.,pt+0.00001};
+ Float_t xval2[2] = {2.+addbin,pt+0.00001};
+
+ TProfile *proTPCpi = pidqa->GetV2(0,xval,xval2);
+ TProfile *proTOFpi = pidqa->GetV2(1,xval,xval2);
+ TProfile *proTPCka = pidqa->GetV2(2,xval,xval2);
+ TProfile *proTOFka = pidqa->GetV2(3,xval,xval2);
+ TProfile *proTPCpr = pidqa->GetV2(4,xval,xval2);
+ TProfile *proTOFpr = pidqa->GetV2(5,xval,xval2);
+
+ proTPCpi->SetName("hPiTPC");
+ proTOFpi->SetName("hPiTOF");
+ proTPCka->SetName("hKaTPC");
+ proTOFka->SetName("hKaTOF");
+ proTPCpr->SetName("hPrTPC");
+ proTOFpr->SetName("hPrTOF");
+
+ proTPCpi->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{#pi} (a.u)");
+ proTOFpi->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{#pi} (ps)");
+ proTPCka->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{K} (a.u)");
+ proTOFka->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{K} (ps)");
+ proTPCpr->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{p} (a.u)");
+ proTOFpr->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{p} (ps)");
+
+ TCanvas *c = new TCanvas("cVproj","cVproj");
+ c->Divide(2,3);
+ c->cd(1);
+ proTPCpi->Draw();
+ c->cd(2);
+ proTOFpi->Draw();
+ c->cd(3);
+ proTPCka->Draw();
+ c->cd(4);
+ proTOFka->Draw();
+ c->cd(5);
+ proTPCpr->Draw();
+ c->cd(6);
+ proTOFpr->Draw();
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff