ClassImp(AliAnalysisTaskRhoVnModulation)
AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation() : AliAnalysisTaskEmcalJet("AliAnalysisTaskRhoVnModulation", kTRUE),
- fDebug(0), fLocalInit(0), fAttachToEvent(kTRUE), fSemiCentralInclusive(kFALSE), fFillHistograms(kTRUE), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(-1.), fReduceBinsYByFactor(-1.), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fTracksCont(0), fJetsCont(0), fUseScaledRho(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fUsePtWeightErrorPropagation(kTRUE), fDetectorType(kTPC), fFitModulationOptions("QWLI"), fRunModeType(kGrid), fDataType(kESD), fCollisionType(kPbPb), fRandom(0), fRunNumber(-1), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fFitControl(0), fMinPvalue(0.01), fMaxPvalue(1), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.), fSemiGoodJetMinPhi(0.), fSemiGoodJetMaxPhi(4.), fSemiGoodTrackMinPhi(0.), fSemiGoodTrackMaxPhi(4.), fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvaluePDF(0), fHistPvalueCDF(0), fHistRhoStatusCent(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fMaxCones(-1), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kTRUE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), fSubtractJetPt(kFALSE), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfV2(0), fProfV2Cumulant(0), fProfV3(0), fProfV3Cumulant(0), fHistPsiControl(0), fHistPsiSpread(0), fHistPsiVZEROA(0), fHistPsiVZEROC(0), fHistPsiVZERO(0), fHistPsiTPC(0), fHistPsiVZEROAV0M(0), fHistPsiVZEROCV0M(0), fHistPsiVZEROVV0M(0), fHistPsiTPCiV0M(0), fHistPsiVZEROATRK(0), fHistPsiVZEROCTRK(0), fHistPsiVZEROTRK(0), fHistPsiTPCTRK(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
+ fDebug(0), fRunToyMC(kFALSE), fLocalInit(0), fAttachToEvent(kTRUE), fSemiCentralInclusive(kFALSE), fFillHistograms(kTRUE), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(-1.), fReduceBinsYByFactor(-1.), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fTracksCont(0), fJetsCont(0), fUseScaledRho(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fFitGoodnessTest(kKolmogorov), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fUsePtWeightErrorPropagation(kTRUE), fDetectorType(kTPC), fFitModulationOptions("QWLI"), fRunModeType(kGrid), fDataType(kESD), fCollisionType(kPbPb), fRandom(0), fRunNumber(-1), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fFitControl(0), fMinPvalue(0.01), fMaxPvalue(1), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.), fSemiGoodJetMinPhi(0.), fSemiGoodJetMaxPhi(4.), fSemiGoodTrackMinPhi(0.), fSemiGoodTrackMaxPhi(4.), fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvalueCDFROOT(0), fHistPvalueCDF(0), fHistKolmogorovTest(0), fHistRhoStatusCent(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fMaxCones(-1), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kTRUE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), fSubtractJetPt(kFALSE), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfV2(0), fProfV2Cumulant(0), fProfV3(0), fProfV3Cumulant(0), fHistPsiControl(0), fHistPsiSpread(0), fHistPsiVZEROA(0), fHistPsiVZEROC(0), fHistPsiVZERO(0), fHistPsiTPC(0), fHistPsiVZEROAV0M(0), fHistPsiVZEROCV0M(0), fHistPsiVZEROVV0M(0), fHistPsiTPCiV0M(0), fHistPsiVZEROATRK(0), fHistPsiVZEROCTRK(0), fHistPsiVZEROTRK(0), fHistPsiTPCTRK(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
for(Int_t i(0); i < 10; i++) {
fProfV2Resolution[i] = 0;
fProfV3Resolution[i] = 0;
}
//_____________________________________________________________________________
AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation(const char* name, runModeType type) : AliAnalysisTaskEmcalJet(name, kTRUE),
- fDebug(0), fLocalInit(0), fAttachToEvent(kTRUE), fSemiCentralInclusive(kFALSE), fFillHistograms(kTRUE), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(-1.), fReduceBinsYByFactor(-1.), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fTracksCont(0), fJetsCont(0), fUseScaledRho(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fUsePtWeightErrorPropagation(kTRUE), fDetectorType(kTPC), fFitModulationOptions("QWLI"), fRunModeType(type), fDataType(kESD), fCollisionType(kPbPb), fRandom(0), fRunNumber(-1), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fFitControl(0), fMinPvalue(0.01), fMaxPvalue(1), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.), fSemiGoodJetMinPhi(0.), fSemiGoodJetMaxPhi(4.), fSemiGoodTrackMinPhi(0.), fSemiGoodTrackMaxPhi(4.), fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvaluePDF(0), fHistPvalueCDF(0), fHistRhoStatusCent(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fMaxCones(-1), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kTRUE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), fSubtractJetPt(kFALSE), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfV2(0), fProfV2Cumulant(0), fProfV3(0), fProfV3Cumulant(0), fHistPsiControl(0), fHistPsiSpread(0), fHistPsiVZEROA(0), fHistPsiVZEROC(0), fHistPsiVZERO(0), fHistPsiTPC(0), fHistPsiVZEROAV0M(0), fHistPsiVZEROCV0M(0), fHistPsiVZEROVV0M(0), fHistPsiTPCiV0M(0), fHistPsiVZEROATRK(0), fHistPsiVZEROCTRK(0), fHistPsiVZEROTRK(0), fHistPsiTPCTRK(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
+ fDebug(0), fRunToyMC(kFALSE), fLocalInit(0), fAttachToEvent(kTRUE), fSemiCentralInclusive(kFALSE), fFillHistograms(kTRUE), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(-1.), fReduceBinsYByFactor(-1.), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fTracksCont(0), fJetsCont(0), fUseScaledRho(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fFitGoodnessTest(kKolmogorov), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fUsePtWeightErrorPropagation(kTRUE), fDetectorType(kTPC), fFitModulationOptions("QWLI"), fRunModeType(type), fDataType(kESD), fCollisionType(kPbPb), fRandom(0), fRunNumber(-1), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fFitControl(0), fMinPvalue(0.01), fMaxPvalue(1), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.), fSemiGoodJetMinPhi(0.), fSemiGoodJetMaxPhi(4.), fSemiGoodTrackMinPhi(0.), fSemiGoodTrackMaxPhi(4.), fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvalueCDFROOT(0), fHistPvalueCDF(0), fHistKolmogorovTest(0), fHistRhoStatusCent(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fMaxCones(-1), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kTRUE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), fSubtractJetPt(kFALSE), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfV2(0), fProfV2Cumulant(0), fProfV3(0), fProfV3Cumulant(0), fHistPsiControl(0), fHistPsiSpread(0), fHistPsiVZEROA(0), fHistPsiVZEROC(0), fHistPsiVZERO(0), fHistPsiTPC(0), fHistPsiVZEROAV0M(0), fHistPsiVZEROCV0M(0), fHistPsiVZEROVV0M(0), fHistPsiTPCiV0M(0), fHistPsiVZEROATRK(0), fHistPsiVZEROCTRK(0), fHistPsiVZEROTRK(0), fHistPsiTPCTRK(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
for(Int_t i(0); i < 10; i++) {
fProfV2Resolution[i] = 0;
fProfV3Resolution[i] = 0;
fProfV3Resolution[i]->GetXaxis()->SetBinLabel(11, "<cos(3(#Psi_{TPC_A} - #Psi_{TPC_B}))>");
fOutputList->Add(fProfV3Resolution[i]);
}
- // cdf and pdf of chisquare distribution
- fHistPvaluePDF = BookTH1F("fHistPvaluePDF", "PDF #chi^{2}", 500, 0, 1);
- fHistPvalueCDF = BookTH1F("fHistPvalueCDF", "CDF #chi^{2}", 500, 0, 1);
- fHistRhoStatusCent = BookTH2F("fHistRhoStatusCent", "centrality", "status [-1=lin was better, 0=ok, 1 = failed]", 101, -1, 100, 3, -1.5, 1.5);
- // vn profile
+ // vn profile
Float_t temp[fCentralityClasses->GetSize()];
for(Int_t i(0); i < fCentralityClasses->GetSize(); i++) temp[i] = fCentralityClasses->At(i);
fProfV2 = new TProfile("fProfV2", "fProfV2", fCentralityClasses->GetSize()-1, temp);
if(fUserSuppliedR3) fOutputList->Add(fUserSuppliedR3);
// increase readability of output list
fOutputList->Sort();
+ // cdf and pdf of chisquare distribution
+ fHistPvalueCDF = BookTH1F("fHistPvalueCDF", "CDF #chi^{2}", 50, 0, 1);
+ fHistKolmogorovTest = BookTH1F("fHistKolmogorovTest", "KolmogorovTest", 50, 0, 1);
+ fHistPvalueCDFROOT = BookTH1F("fHistPvalueCDFROOT", "CDF #chi^{2} ROOT", 50, 0, 1);
+ fHistRhoStatusCent = BookTH2F("fHistRhoStatusCent", "centrality", "status [-1=lin was better, 0=ok, 1 = failed]", 101, -1, 100, 3, -1.5, 1.5);
+
PostData(1, fOutputList);
switch (fRunModeType) {
} break;
default : break;
}
+ if(fRunToyMC) {
+ // toy mc, just here to check procedure, azimuthal profile is filled from hypothesis so p-value distribution should be flat
+ Int_t _bins = _tempSwap.GetXaxis()->GetNbins();
+ TF1* _tempFit = new TF1("temp_fit_kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi());
+ _tempFit->SetParameter(0, fFitModulation->GetParameter(0)); // normalization
+ _tempFit->SetParameter(3, 0.1); // v2
+ _tempFit->FixParameter(1, 1.); // constant
+ _tempFit->FixParameter(2, 2.); // constant
+ _tempFit->FixParameter(5, 3.); // constant
+ _tempFit->FixParameter(4, fFitModulation->GetParameter(4));
+ _tempFit->FixParameter(6, fFitModulation->GetParameter(6));
+ _tempFit->SetParameter(7, 0.1); // v3
+ _tempSwap.Reset(); // rese bin content
+ for(int _binsI = 0; _binsI < _bins*_bins; _binsI++) _tempSwap.Fill(_tempFit->GetRandom());
+ }
_tempSwap.Fit(fFitModulation, fFitModulationOptions.Data(), "", lowBound, upBound);
// the quality of the fit is evaluated from 1 - the cdf of the chi square distribution
- Double_t CDF(1.-ChiSquareCDF(fFitModulation->GetNDF(), fFitModulation->GetChisquare()));
+ // three methods are available, all with their drawbacks. all are stored, one is selected to do the cut
+ Double_t CDF(1.-ChiSquareCDF(fFitModulation->GetNDF(), ChiSquare(_tempSwap, fFitModulation)));
+ Double_t CDFROOT(1.-ChiSquareCDF(fFitModulation->GetNDF(), fFitModulation->GetChisquare()));
+ Double_t CDFKolmogorov(KolmogorovTest(_tempSwap, fFitModulation));
fHistPvalueCDF->Fill(CDF);
+ fHistPvalueCDFROOT->Fill(CDFROOT);
+ fHistKolmogorovTest->Fill(CDFKolmogorov);
+ // variable CDF is used for making cuts, so we fill it with the selected p-value
+ switch (fFitGoodnessTest) {
+ case kChi2ROOT : {
+ CDF = CDFROOT;
+ } break;
+ case kChi2Poisson : break; // CDF is already CDF
+ case kKolmogorov : {
+ CDF = CDFKolmogorov;
+ } break;
+ default: break;
+ }
+
if(fFitControl) {
// as an additional quality check, see if fitting a control fit has a higher significance
_tempSwap.Fit(fFitControl, fFitModulationOptions.Data(), "", lowBound, upBound);
- Double_t CDFControl(1.-ChiSquareCDF(fFitControl->GetNDF(), fFitControl->GetChisquare()));
+ Double_t CDFControl(-1.);
+ switch (fFitGoodnessTest) {
+ case kChi2ROOT : {
+ CDFControl = 1.-ChiSquareCDF(fFitControl->GetNDF(), fFitModulation->GetChisquare());
+ } break;
+ case kChi2Poisson : {
+ CDFControl = 1.-ChiSquareCDF(fFitControl->GetNDF(), ChiSquare(_tempSwap, fFitModulation));
+ } break;
+ case kKolmogorov : {
+ CDFControl = KolmogorovTest(_tempSwap, fFitControl);
+ } break;
+ default: break;
+ }
if(CDFControl > CDF) {
CDF = -1.; // control fit is more significant, so throw out the 'old' fit
fHistRhoStatusCent->Fill(fCent, -1);
public:
// enumerators
enum fitModulationType { kNoFit, kV2, kV3, kCombined, kFourierSeries, kIntegratedFlow, kQC2, kQC4 }; // fit type
+ enum fitGoodnessTest { kChi2ROOT, kChi2Poisson, kKolmogorov, kKolmogorovTOY, kLinearFit };
enum collisionType { kPbPb, kPythia }; // collision type
enum qcRecovery { kFixedRho, kNegativeVn, kTryFit }; // how to deal with negative cn value for qcn value
enum runModeType { kLocal, kGrid }; // run mode type
return -999; }
// note that the cdf of the chisquare distribution is the normalized lower incomplete gamma function
/* inline */ Double_t ChiSquareCDF(Int_t ndf, Double_t x) const { return TMath::Gamma(ndf/2., x/2.); }
+ /* inline */ Double_t ChiSquare(TH1& histo, TF1* func) const {
+ // evaluate the chi2 using a poissonian error estimate on bins
+ Double_t chi2(0.);
+ for(Int_t i(0); i < histo.GetXaxis()->GetNbins(); i++) {
+ if(histo.GetBinContent(i+1) <= 0.) continue;
+ chi2 += TMath::Power((histo.GetBinContent(i+1)-func->Eval(histo.GetXaxis()->GetBinCenter(1+i))), 2)/histo.GetBinContent(i+1);
+ }
+ return chi2;
+ }
+ /* inline*/ Double_t KolmogorovTest(TH1F& histo, TF1* func) const {
+ // return the probability from a Kolmogorov test
+ 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");
+ return histo.TH1::KolmogorovTest((&test));
+ }
+
// setters - analysis setup
void SetDebugMode(Int_t d) {fDebug = d;}
+ void SetRunToyMC(Bool_t t) {fRunToyMC = t; }
void SetAttachToEvent(Bool_t b) {fAttachToEvent = b;}
void SetSemiCentralInclusive(Bool_t b) {fSemiCentralInclusive = b;}
void SetFillHistograms(Bool_t b) {fFillHistograms = b;}
void SetUseControlFit(Bool_t c);
void SetModulationFitMinMaxP(Float_t m, Float_t n) {fMinPvalue = m; fMaxPvalue = n; }
void SetModulationFitType(fitModulationType type) {fFitModulationType = type; }
+ void SetGoodnessTest(fitGoodnessTest test) {fFitGoodnessTest = test; }
void SetQCnRecoveryType(qcRecovery type) {fQCRecovery = type; }
void SetModulationFitOptions(TString opt) {fFitModulationOptions = opt; }
void SetReferenceDetector(detectorType type) {fDetectorType = type; }
void SetCollisionType(collisionType type) {fCollisionType = type; }
- void SetUsePtWeight(Bool_t w) {fUsePtWeight = w; }
+ void SetUsePtWeight(Bool_t w) {
+ fUsePtWeight = w;
+ if(!fUsePtWeight) fUsePtWeightErrorPropagation = kFALSE; }
void SetUsePtWeightErrorPropagation(Bool_t w) {fUsePtWeightErrorPropagation = w; }
void SetRunModeType(runModeType type) {fRunModeType = type; }
void SetAbsVertexZ(Float_t v) {fAbsVertexZ = v; }
private:
// analysis flags and settings
Int_t fDebug; // debug level (0 none, 1 fcn calls, 2 verbose)
+ Bool_t fRunToyMC; // run toy mc for fit routine
Bool_t fLocalInit; //! is the analysis initialized?
Bool_t fAttachToEvent; // attach local rho to the event
Bool_t fSemiCentralInclusive; // semi central inclusive event selection
Int_t fNAcceptedTracks; //! number of accepted tracks
Int_t fNAcceptedTracksQCn; //! accepted tracks for QCn
fitModulationType fFitModulationType; // fit modulation type
+ fitGoodnessTest fFitGoodnessTest; // fit goodness test type
qcRecovery fQCRecovery; // recovery type for e-by-e qc method
Bool_t fUsePtWeight; // use dptdphi instead of dndphi
Bool_t fUsePtWeightErrorPropagation; // recalculate the bin errors in case of pt weighting
TH1F* fHistVertexz; //! accepted verte
TH2F* fHistRunnumbersPhi; //! run numbers averaged phi
TH2F* fHistRunnumbersEta; //! run numbers averaged eta
- TH1F* fHistPvaluePDF; //! pdf value of chisquare p
+ TH1F* fHistPvalueCDFROOT; //! pdf value of chisquare p
TH1F* fHistPvalueCDF; //! cdf value of chisquare p
+ TH1F* fHistKolmogorovTest; //! KolmogorovTest value
TH2F* fHistRhoStatusCent; //! status of rho as function of centrality
// general settings
Float_t fMinDisanceRCtoLJ; // min distance between rc and leading jet