fResponseMaker (new AliAnaChargedJetResponseMaker()),
fRMS (kTRUE),
fSymmRMS (kTRUE),
+ fRho0 (kFALSE),
fPower (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
fSaveFull (kTRUE),
fActiveString (""),
}
// extract the spectra
TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
+ if(fRho0) spectrumName = Form("fHistJetPsi2PtRho0_%i", fCentralityArray->At(0));
if(!fInputList->FindObject(spectrumName.Data())) {
printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
return kFALSE;
}
// extract the delta pt matrices
TString deltaptName("");
- deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
+ if(!fRho0) {
+ deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
+ } else {
+ deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJRho0_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2Rho0_%i", fCentralityArray->At(0));
+ }
fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
if(!fDeltaPtDeltaPhi) {
printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
// clear minuit state to avoid constraining the fit with the results of the previous iteration
for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
}
- // extract the spectra
+ // extract the spectra
TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
+ if(fRho0) spectrumName = Form("fHistJetPsi2PtRho0_%i", fCentralityArray->At(0));
fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
if(!fJetPtDeltaPhi) {
printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), low, up, "e");
// extract the delta pt matrices
TString deltaptName("");
- deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
+ if(!fRho0) {
+ deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
+ } else {
+ deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJRho0_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2Rho0_%i", fCentralityArray->At(0));
+ }
fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
if(!fDeltaPtDeltaPhi) {
printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
h->SetMarkerStyle(8);
h->SetMarkerSize(1);
} break;
+ case kDeltaPhi : {
+ h->GetYaxis()->SetTitle("[counts]");
+ h->GetXaxis()->SetTitle("#Delta #phi");
+ }
default : break;
}
}
Int_t low[] = {1, 6, 11, 16, 21, 26, 31, 36};
Int_t up[] = {5, 10, 15, 20, 25, 30, 35, 40};
TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
- TH1D* dPtdPhi[8];
- for(Int_t i(0); i < 8; i++) dPtdPhi[i] = new TH1D(Form("dPtdPhi_%i", i), Form("dPtdPhi_%i", i), 8, 0, TMath::Pi());
+ const Int_t ptBins(fBinsTrue->GetSize()-1);
+ const Int_t dPhiBins(8);
+ TH1D* dPtdPhi[fBinsTrue->GetSize()];
+ for(Int_t i(0); i < ptBins; i++) dPtdPhi[i] = new TH1D(Form("dPtdPhi_%i", i), Form("dPtdPhi_%i", i), dPhiBins, 0, TMath::Pi());
- for(Int_t i(0); i < 8; i++) {
+ // for the output initialize a canvas
+ TCanvas* v2Fits(new TCanvas("v2 fits", "v2 fits"));
+ v2Fits->Divide(4, TMath::Floor((1+ptBins)/(float)4)+(((1+ptBins)%4)>0));
+
+ for(Int_t i(0); i < dPhiBins; i++) {
// 1) manipulation of input histograms
// check if the input variables are present
if(!PrepareForUnfolding(low[i], up[i])) return;
measuredJetSpectrumTrueBinsIn,
TString("dPtdPhiUnfolding"),
jetFindingEfficiency);
- if(i==5) {
+ // arbitrarily save one of the full outputs (same for all dphi bins, avoid duplicates)
+ if(i+1 == ptBins) {
resizedResponseIn->SetNameTitle(Form("ResponseMatrix_%s", stringArray[i].Data()), Form("response matrix %s", stringArray[i].Data()));
resizedResponseIn->SetXTitle("p_{T, jet}^{true} [GeV/c]");
resizedResponseIn->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
TH1D* dud(ProtectHeap(unfoldedJetSpectrumIn, kTRUE, stringArray[i]));;
Double_t integralError(0);
- for(Int_t j(0); j < 6; j++) {
+ // at this point in the code, the spectrum has been unfolded in a certain region of dPhi space
+ // next step is splitting it in pt space as well to estimate the yield differentially in pt
+ for(Int_t j(0); j < ptBins; j++) {
// get the integrated
- Double_t integral(dud->IntegralAndError(2*j+1, 2*j+3, integralError));
+ Double_t integral(dud->IntegralAndError(j+1, j+2, integralError));
dPtdPhi[j]->SetBinContent(i+1, integral);
dPtdPhi[j]->SetBinError(i+1, integralError);
}
}
TF1* fourier = new TF1("fourier", "[0]*(1.+0.5*[1]*(TMath::Cos(2.*x)))", 0, TMath::Pi());
TH1D* v2(new TH1D("v2FromFit", "v2FromFit", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
- for(Int_t i(0); i < 6; i++) {
+ for(Int_t i(0); i < ptBins; i++) {
+ v2Fits->cd(i+1);
dPtdPhi[i]->Fit(fourier, "VI");
+ Style(gPad, "PEARSON");
+ Style(dPtdPhi[i], kBlue, kDeltaPhi);
+ dPtdPhi[i]->DrawCopy();
+ AliJetFlowTools::AddText(
+ TString(Form("%.2f #LT p_{T} #LT %.2f", fBinsTrue->At(i), fBinsTrue->At(i+1))),
+ kBlack,
+ .38,
+ .56,
+ .62,
+ .65
+ );
v2->SetBinContent(1+i, fourier->GetParameter(1));
v2->SetBinError(1+i, fourier->GetParError(1));
- dPtdPhi[i]->Write();
}
- v2->Write();
+ v2Fits->cd(1+ptBins);
+ Style(gPad, "PEARSON");
+ Style(v2, kBlack, kV2, kTRUE);
+ v2->DrawCopy();
+ v2Fits->Write();
}
//_____________________________________________________________________________
void AliJetFlowTools::ReplaceBins(TArrayI* array, TGraphErrors* graph) {
void AliAnalysisTaskJetV2::CalculateEventPlaneVZERO(Double_t vzero[2][2]) const
{
// get the vzero event plane
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
if(fUseV0EventPlaneFromHeader) { // use the vzero from the header
Double_t a(0), b(0), c(0), d(0), e(0), f(0), g(0), h(0);
vzero[0][0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, a, b);
void AliAnalysisTaskJetV2::CalculateEventPlaneCombinedVZERO(Double_t* comb) const
{
// grab the combined vzero event plane
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
Double_t a(0), b(0), c(0), d(0);
comb[0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, a, b);
comb[1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, c, d);
if(GetParticleContainer()->GetParticlePhiMin() > lowBound) lowBound = GetParticleContainer()->GetParticlePhiMin();
if(GetParticleContainer()->GetParticlePhiMax() < upBound) upBound = GetParticleContainer()->GetParticlePhiMax();
- fHistSwap->Reset(); // clear the histogram
+ fHistSwap->Reset(); // clear the histogram
TH1F _tempSwap; // on stack for quick access
TH1F _tempSwapN; // on stack for quick access, bookkeeping histogram
if(fRebinSwapHistoOnTheFly) {
// enumerators
enum fitModulationType { kNoFit, kV2, kV3, kCombined, kFourierSeries, kIntegratedFlow, kQC2, kQC4 }; // fit type
enum fitGoodnessTest { kChi2ROOT, kChi2Poisson, kKolmogorov, kKolmogorovTOY, kLinearFit };
- enum collisionType { kPbPb, kPythia, kPbPb10h }; // collision type
+ enum collisionType { kPbPb, kPythia, kPbPb10h, kPbPb11h }; // collision type, kPbPb = 11h, kept for backward compatibilitiy
enum qcRecovery { kFixedRho, kNegativeVn, kTryFit }; // how to deal with negative cn value for qcn value
enum runModeType { kLocal, kGrid }; // run mode type
enum dataType { kESD, kAOD, kESDMC, kAODMC }; // data type
- enum detectorType { kTPC, kVZEROA, kVZEROC, kVZEROComb}; // detector that was used
+ enum detectorType { kTPC, kVZEROA, kVZEROC, kVZEROComb}; // detector that was used for event plane
enum analysisType { kCharged, kFull }; // analysis type
// constructors, destructor
AliAnalysisTaskJetV2();
}
/* inline*/ Double_t KolmogorovTest(TH1F& histo, TF1* func) const {
// return the probability from a Kolmogorov test
+ return .5;
TH1F test(histo); // stack copy of test statistic
for(Int_t i(0); i < test.GetXaxis()->GetNbins(); i++) test.SetBinContent(i+1, func->Eval(test.GetXaxis()->GetBinCenter(1+i)));
if(fFitGoodnessTest == kKolmogorovTOY) return histo.TH1::KolmogorovTest((&test), "X");