from Redmer

diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.cxx
index b681f931b4a8c194103a1b0a51c496618adaa796..2cc31121e5f776175b28702630a6862cf387f9e8 100644 (file)
--- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.cxx
+++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.cxx
@@ -61,7 +61,7 @@ using namespace std;
 ClassImp(AliAnalysisTaskRhoVnModulation)
 
 AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation() : AliAnalysisTaskEmcalJet("AliAnalysisTaskRhoVnModulation", kTRUE), 
-    fDebug(0), fInitialized(0), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(1), fReduceBinsYByFactor(1), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fUsePtWeight(kTRUE), fDetectorType(kTPC), fFitModulationOptions("Q"), fRunModeType(kGrid), fDataType(kESD), fRandom(0), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fMinPvalue(0), fMaxPvalue(1), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.), fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvaluePDF(0), fHistPvalueCDF(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kFALSE), fSetPtSub(kFALSE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), 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), fHistPsiTPC(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
+    fDebug(0), fInitialized(0), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(1), fReduceBinsYByFactor(1), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fDetectorType(kTPC), fFitModulationOptions("Q"), fRunModeType(kGrid), fDataType(kESD), fRandom(0), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(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.), fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvaluePDF(0), fHistPvalueCDF(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kFALSE), fSetPtSub(kFALSE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), 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), fHistPsiTPC(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
     for(Int_t i(0); i < 10; i++) {
         fProfV2Resolution[i] = 0;
         fProfV3Resolution[i] = 0;
@@ -120,11 +120,12 @@ AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation() : AliAnalysisTa
 }
 //_____________________________________________________________________________
 AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation(const char* name, runModeType type) : AliAnalysisTaskEmcalJet(name, kTRUE),
-  fDebug(0), fInitialized(0), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(1), fReduceBinsYByFactor(1), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fUsePtWeight(kTRUE), fDetectorType(kTPC), fFitModulationOptions("Q"), fRunModeType(type), fDataType(kESD), fRandom(0), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fMinPvalue(0), fMaxPvalue(1), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.),  fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvaluePDF(0), fHistPvalueCDF(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kFALSE), fSetPtSub(kFALSE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), 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), fHistPsiTPC(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
+  fDebug(0), fInitialized(0), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(1), fReduceBinsYByFactor(1), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fDetectorType(kTPC), fFitModulationOptions("Q"), fRunModeType(type), fDataType(kESD), fRandom(0), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(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.),  fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvaluePDF(0), fHistPvalueCDF(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kFALSE), fSetPtSub(kFALSE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), 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), fHistPsiTPC(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
     for(Int_t i(0); i < 10; i++) {
         fProfV2Resolution[i] = 0;
         fProfV3Resolution[i] = 0;
         fHistPicoTrackPt[i] = 0;
+        fHistPicoTrackMult[i] = 0;
         fHistPicoCat1[i] = 0;
         fHistPicoCat2[i] = 0;
         fHistPicoCat3[i] = 0;
@@ -299,6 +300,7 @@ void AliAnalysisTaskRhoVnModulation::UserCreateOutputObjects()
     // pico track kinematics
     for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) { 
         fHistPicoTrackPt[i] =          BookTH1F("fHistPicoTrackPt", "p_{t} [GeV/c]", 100, 0, 50, i);
+        fHistPicoTrackMult[i] =        BookTH1F("fHistPicoTrackMult", "multiplicity", 100, 0, 5000, i);
         if(fFillQAHistograms) {
             fHistPicoCat1[i] =             BookTH2F("fHistPicoCat1", "#eta", "#phi", 50, -1, 1, 50, 0, TMath::TwoPi(), i);
             fHistPicoCat2[i] =             BookTH2F("fHistPicoCat2", "#eta", "#phi", 50, -1, 1, 50, 0, TMath::TwoPi(), i);
@@ -497,7 +499,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
     }
     switch (fFitModulationType) { // do the fits
         case kNoFit : { fFitModulation->FixParameter(0, RhoVal()); } break;
-        case kV2 : {
+        case kV2 : {    // only v2
             if(CorrectRho(psi2, psi3)) {
                 fProfV2->Fill(fCent, fFitModulation->GetParameter(3));
                 if(fUserSuppliedR2) {
@@ -507,7 +509,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
                 CalculateEventPlaneResolution(vzero, tpc);
             }
         } break;
-        case kV3 : {
+        case kV3 : {    // only v3
             if(CorrectRho(psi2, psi3)) {
                 if(fUserSuppliedR3) {
                     Double_t r(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
@@ -517,54 +519,42 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
                 CalculateEventPlaneResolution(vzero, tpc);
             }
         } break;
-        case kQC2 : {
-            Bool_t QC2(CorrectRho(psi2, psi3));
-            if(fUserSuppliedR2 && fUserSuppliedR3) {
-                // note for the qc method, resolution is REVERSED to go back to v2obs
-                Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
-                Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
-                if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
-                if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(3)*r3);
-            }
-            if (fUsePtWeight) { // use weighted weights
-                fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3), QCnM11());
-                fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7), QCnM11());
-                if(QC2) {        // how to deal with negative results from the cumulants ?
-                    fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.5), QCnM11());
-                    fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.5), QCnM11()); 
+        case kQC2 : {   // qc2 analysis
+            if(CorrectRho(psi2, psi3)) {
+                if(fUserSuppliedR2 && fUserSuppliedR3) {
+                    // note for the qc method, resolution is REVERSED to go back to v2obs
+                    Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
+                    Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
+                    if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
+                    if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(3)*r3);
                 }
-            } else {
-                fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3), QCnM()*(QCnM()-1));
-                fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7), QCnM()*(QCnM()-1));
-                if(QC2) {
-                    fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.5), QCnM()*(QCnM()-1));
-                    fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.5), QCnM()*(QCnM()-1));
+                if (fUsePtWeight) { // use weighted weights
+                    Double_t dQCnM11 = (fNoEventWeightsForQC) ? 1. : QCnM11();
+                    fProfV2->Fill(fCent, fFitModulation->GetParameter(3), dQCnM11);
+                    fProfV3->Fill(fCent, fFitModulation->GetParameter(7), dQCnM11); 
+                } else {
+                    Double_t dQCnM = (fNoEventWeightsForQC) ? 2. : QCnM();
+                    fProfV2->Fill(fCent, fFitModulation->GetParameter(3), dQCnM*(dQCnM-1));
+                    fProfV3->Fill(fCent, fFitModulation->GetParameter(7), dQCnM*(dQCnM-1));
                 }
+                CalculateEventPlaneResolution(vzero, tpc);
             }
-            CalculateEventPlaneResolution(vzero, tpc);
         } break;
         case kQC4 : {
-            Bool_t QC4(CorrectRho(psi2, psi3));
-            if(fUserSuppliedR2 && fUserSuppliedR3) {
-                // note for the qc method, resolution is REVERSED to go back to v2obs
-                Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
-                Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
-                if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
-                if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(3)*r3);
-            }
-            if (fUsePtWeight) { // use weighted weights
-                fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3)/*, QCnM1111()*/);
-                fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7)/*, QCnM1111()*/);
-                if(QC4) {        // how to deal with negative values of cumulants ?
+            if(CorrectRho(psi2, psi3)) {
+                if(fUserSuppliedR2 && fUserSuppliedR3) {
+                    // note for the qc method, resolution is REVERSED to go back to v2obs   
+                    Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
+                    Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
+                    if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
+                    if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(3)*r3);
+                }
+                if (fUsePtWeight) { // use weighted weights
                     fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.5)/*, QCnM1111()*/);
                     fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.5)/*, QCnM1111()*/); 
-                }
-            } else {
-                fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
-                fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7)/*,  QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
-                if(QC4) {
-                    fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.25)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
-                    fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.25)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
+                } else {
+                    fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.5)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
+                    fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.5)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
                 }
             }
             CalculateEventPlaneResolution(vzero, tpc);
@@ -878,6 +868,46 @@ Double_t AliAnalysisTaskRhoVnModulation::QCnM1111() {
     return (QCnS(4,1)-6*QCnS(1,2)*QCnS(2,1)+8*QCnS(1,3)*QCnS(1,1)+3*QCnS(2,2)-6*QCnS(1,4));
 }
 //_____________________________________________________________________________
+Bool_t AliAnalysisTaskRhoVnModulation::QCnRecovery(Double_t psi2, Double_t psi3) {
+    // decides how to deal with the situation where c2 or c3 is negative 
+    // returns kTRUE depending on whether or not a modulated rho is used for the jet background
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    if(TMath::AreEqualAbs(fFitModulation->GetParameter(3), .0, 1e-10) && TMath::AreEqualAbs(fFitModulation->GetParameter(7), .0,1e-10)) {
+        fFitModulation->SetParameter(7, 0);
+        fFitModulation->SetParameter(3, 0);
+        fFitModulation->SetParameter(0, RhoVal());
+        return kTRUE;   // v2 and v3 have physical null values
+    }
+    switch (fQCRecovery) {
+        case kFixedRho : {      // roll back to the original rho
+           fFitModulation->SetParameter(7, 0);
+           fFitModulation->SetParameter(3, 0);
+           fFitModulation->SetParameter(0, RhoVal());
+           return kFALSE;       // rho is forced to be fixed
+        }
+        case kNegativeVn : {
+           Double_t c2(fFitModulation->GetParameter(3));
+           Double_t c3(fFitModulation->GetParameter(7));
+           if( c2 < 0 ) c2 = -1.*TMath::Sqrt(-1.*c2);
+           if( c3 < 0 ) c3 = -1.*TMath::Sqrt(-1.*c3);
+           fFitModulation->SetParameter(3, c2);
+           fFitModulation->SetParameter(7, c3);
+           return kTRUE;        // is this a physical quantity ?
+        }
+        case kTryFit : {
+           fitModulationType tempType(fFitModulationType);  // store temporarily
+           fFitModulationType = kCombined;
+           fFitModulation->SetParameter(7, 0);
+           fFitModulation->SetParameter(3, 0);
+           Bool_t pass(CorrectRho(psi2, psi3));         // do the fit and all quality checks
+           fFitModulationType = tempType;               // roll back for next event
+           return pass;
+        }
+        default : return kFALSE;
+    }
+    return kFALSE;
+}
+//_____________________________________________________________________________
 Bool_t AliAnalysisTaskRhoVnModulation::CorrectRho(Double_t psi2, Double_t psi3) 
 {
     // get rho' -> rho(phi)
@@ -886,6 +916,9 @@ Bool_t AliAnalysisTaskRhoVnModulation::CorrectRho(Double_t psi2, Double_t psi3)
     //      in case of cumulants, both cumulants and vn values are stored. in both cases, v2 and v3
     //      are expected. a check is performed to see if rho has no negative local minimum
     //      for full description, see Phys. Rev. C 83, 044913
+    //      since the cn distribution has negative values, vn = sqrt(cn) can be imaginary sometimes
+    //      in this case one can either roll back to the 'original' rixed rho, do a fit for vn or take use
+    //      vn = - sqrt(|cn|) 
     //  [2] fitting a fourier expansion to the de/dphi distribution
     //      the fit can be done with either v2, v3 or a combination.
     //      in all cases, a cut can be made on the p-value of the chi-squared value of the fit
@@ -895,28 +928,55 @@ Bool_t AliAnalysisTaskRhoVnModulation::CorrectRho(Double_t psi2, Double_t psi3)
         case kQC2 : {
             fFitModulation->FixParameter(4, psi2); 
             fFitModulation->FixParameter(6, psi3);
-            fFitModulation->FixParameter(3, CalculateQC2(2));
+            fFitModulation->FixParameter(3, CalculateQC2(2));   // set here with cn, vn = sqrt(cn)
             fFitModulation->FixParameter(7, CalculateQC2(3));
-            if(fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) { 
+            // first fill the histos of the raw cumulant distribution
+            if (fUsePtWeight) { // use weighted weights
+                Double_t dQCnM11 = (fNoEventWeightsForQC) ? 1. : QCnM11();
+                fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3), dQCnM11);
+                fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7), dQCnM11);
+            } else {
+                Double_t dQCnM = (fNoEventWeightsForQC) ? 2. : QCnM();
+                fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3), dQCnM*(dQCnM-1));
+                fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7), dQCnM*(dQCnM-1));
+            }
+            // then see if one of the cn value is larger than zero and vn is readily available
+            if(fFitModulation->GetParameter(3) > 0 && fFitModulation->GetParameter(7) > 0) {
+                fFitModulation->FixParameter(3, TMath::Sqrt(fFitModulation->GetParameter(3)));
+                fFitModulation->FixParameter(7, TMath::Sqrt(fFitModulation->GetParameter(7)));
+            } else if (!QCnRecovery(psi2, psi3)) return kFALSE;  // try to recover the cumulant, this will set v2 and v3
+            if(fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) {  // general check 
                 fFitModulation->SetParameter(7, 0);
                 fFitModulation->SetParameter(3, 0);
                 fFitModulation->SetParameter(0, RhoVal());
                 return kFALSE;
             }
-            return (fFitModulation->GetParameter(3) > 0 && fFitModulation->GetParameter(7) > 0) ? kTRUE : kFALSE;
+            return kTRUE;
         } break;
         case kQC4 : {
             fFitModulation->FixParameter(4, psi2); 
             fFitModulation->FixParameter(6, psi3);
-            fFitModulation->FixParameter(3, CalculateQC4(2));
+            fFitModulation->FixParameter(3, CalculateQC4(2));   // set here with cn, vn = sqrt(cn)
             fFitModulation->FixParameter(7, CalculateQC4(3));
-            if(fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) { 
+            // first fill the histos of the raw cumulant distribution
+            if (fUsePtWeight) { // use weighted weights
+                fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3)/*, QCnM1111()*/);
+                fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7)/*, QCnM1111()*/);
+            } else {
+                fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3)/*, QCnM1111()*/);
+                fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7)/*, QCnM1111()*/);
+            }
+            // then see if one of the cn value is larger than zero and vn is readily available
+            if(fFitModulation->GetParameter(3) > 0 && fFitModulation->GetParameter(7) > 0) {
+                fFitModulation->FixParameter(3, TMath::Sqrt(fFitModulation->GetParameter(3)));
+                fFitModulation->FixParameter(7, TMath::Sqrt(fFitModulation->GetParameter(7)));
+            } else if (!QCnRecovery(psi2, psi3)) return kFALSE;  // try to recover the cumulant, this will set v2 and v3
+            if(fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) {  // general check 
                 fFitModulation->SetParameter(7, 0);
                 fFitModulation->SetParameter(3, 0);
                 fFitModulation->SetParameter(0, RhoVal());
                 return kFALSE;
             }
-            return (fFitModulation->GetParameter(3) > 0 && fFitModulation->GetParameter(7) > 0) ? kTRUE : kFALSE;
         } break;
         case kIntegratedFlow : {
             // use v2 and v3 values from an earlier iteration over the data
@@ -992,7 +1052,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::CorrectRho(Double_t psi2, Double_t psi3)
         case kV3 : { 
             fFitModulation->FixParameter(4, psi3); 
         } break;
-        case kCombined : { 
+        case kCombined : {
             fFitModulation->FixParameter(4, psi2); 
             fFitModulation->FixParameter(6, psi3);
         } break;
@@ -1019,9 +1079,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::CorrectRho(Double_t psi2, Double_t psi3)
     _tempSwap.Fit(fFitModulation, fFitModulationOptions.Data(), "", 0, TMath::TwoPi());
     // 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()));
-//    Double_t PDF(ChiSquarePDF(fFitModulation->GetNDF(), fFitModulation->GetChisquare()));
     fHistPvalueCDF->Fill(CDF);
-//    fHistPvaluePDF->Fill(PDF);
     if(CDF > fMinPvalue && CDF < fMaxPvalue && ( fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) > 0)) { // fit quality
         // for LOCAL didactic purposes, save the  best and the worst fits
         // this routine can produce a lot of output histograms (it's not memory 'safe') and will not work on GRID 
@@ -1164,13 +1222,15 @@ void AliAnalysisTaskRhoVnModulation::FillTrackHistograms() const
 {
     // fill track histograms
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
-    Int_t iTracks(fTracks->GetEntriesFast());
+    Int_t iTracks(fTracks->GetEntriesFast()), iAcceptedTracks(0);
     for(Int_t i(0); i < iTracks; i++) {
         AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
         if(!PassesCuts(track)) continue;
+        iAcceptedTracks++;
         fHistPicoTrackPt[fInCentralitySelection]->Fill(track->Pt());
         if(fFillQAHistograms) FillQAHistograms(track);
     }
+    fHistPicoTrackMult[fInCentralitySelection]->Fill(iAcceptedTracks);
 }
 //_____________________________________________________________________________
 void AliAnalysisTaskRhoVnModulation::FillClusterHistograms() const
@@ -1561,3 +1621,23 @@ TH1F* AliAnalysisTaskRhoVnModulation::CorrectForResolutionInt(TH1F* v, detectorT
     return (v->Divide(v, r)) ? v : 0x0;
 }
 //_____________________________________________________________________________
+TH1F* AliAnalysisTaskRhoVnModulation::GetDifferentialQC(TProfile* refCumulants, TProfile* diffCumlants, TArrayD* ptBins, Int_t h)
+{
+    // get differential QC
+    Double_t r(refCumulants->GetBinContent(h-1)); // v2 reference flow
+    if(r > 0) r = TMath::Sqrt(r);
+    TH1F* qc = new TH1F(Form("QC2v%i", h), Form("QC2v%i", h), ptBins->GetSize()-1, ptBins->GetArray());
+    Double_t a(0), b(0), c(0);  // dummy variables
+    for(Int_t i(0); i < ptBins->GetSize(); i++) {
+        if(r > 0) {
+            a = diffCumlants->GetBinContent(1+i);
+            b = diffCumlants->GetBinError(1+i);
+            c = a/r;
+            qc->SetBinContent(1+i, c);
+            (a <= 0 || b <= 0) ? qc->SetBinError(1+i, b) : qc->SetBinError(1+i, TMath::Sqrt(c*c*b*b/(a*a)));
+        }
+    }
+    return qc;
+}
+
+//_____________________________________________________________________________

diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.h b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.h
index b0a22f23c48ecc8a3d5b726769dd359e0cddb9b5..5800e76296875160395b066badbc58fd27d27470 100644 (file)
--- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.h
+++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.h
@@ -24,6 +24,7 @@ class AliAnalysisTaskRhoVnModulation : public AliAnalysisTaskEmcalJet
     public:
          // enumerators
         enum fitModulationType  { kNoFit, kV2, kV3, kCombined, kFourierSeries, kIntegratedFlow, kQC2, kQC4 }; // fit type
+        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};    // detector that was used
@@ -72,6 +73,7 @@ class AliAnalysisTaskRhoVnModulation : public AliAnalysisTaskEmcalJet
         void                    SetFillQAHistograms(Bool_t qa)                  {fFillQAHistograms = qa;}
         void                    SetReduceBinsXYByFactor(Int_t x, Int_t y)       {fReduceBinsXByFactor = x;
                                                                                  fReduceBinsYByFactor = y;}
+        void                    SetNoEventWeightsForQC(Bool_t e)                {fNoEventWeightsForQC = e;}
         void                    SetCentralityClasses(TArrayI* c)                {fCentralityClasses = c;}
         void                    SetPtBinsHybrids(TArrayD* p)                    {fPtBinsHybrids = p;}
         void                    SetPtBinsJets(TArrayD* p)                       {fPtBinsJets = p;}
@@ -87,6 +89,7 @@ class AliAnalysisTaskRhoVnModulation : public AliAnalysisTaskEmcalJet
                                                                                  fFitModulation = fit; }
         void                    SetModulationFitMinMaxP(Float_t m, Float_t n)   {fMinPvalue = m; fMaxPvalue = n; }
         void                    SetModulationFitType(fitModulationType type)    {fFitModulationType = type; }
+        void                    SetQCnRecoveryType(qcRecovery type)             {fQCRecovery = type; }
         void                    SetModulationFitOptions(TString opt)            {fFitModulationOptions = opt; }
         void                    SetReferenceDetector(detectorType type)         {fDetectorType = type; }
         void                    SetUsePtWeight(Bool_t w)                        {fUsePtWeight = w; }
@@ -111,6 +114,7 @@ class AliAnalysisTaskRhoVnModulation : public AliAnalysisTaskEmcalJet
         TProfile*               GetResolutionParameters(Int_t h, Int_t c) const {return (h==2) ? fProfV2Resolution[c] : fProfV3Resolution[c];}
         TList*                  GetOutputList() const                           {return fOutputList;}
         void                    ExecMe()                                        {ExecOnce();}
+        AliAnalysisTaskRhoVnModulation* ReturnMe()                              {return this;}
         // local cuts
         void                    SetLocalJetMinMaxEta(Float_t min, Float_t max)  {fLocalJetMinEta = min; fLocalJetMaxEta = max;}
         void                    SetLocalJetMinMaxEta(Float_t R)                 {fLocalJetMinEta = - 0.9 + R; fLocalJetMaxEta = 0.9 - R; }
@@ -132,6 +136,7 @@ class AliAnalysisTaskRhoVnModulation : public AliAnalysisTaskEmcalJet
         Double_t                QCnM();
         Double_t                QCnM11();
         Double_t                QCnM1111();
+        Bool_t                  QCnRecovery(Double_t psi2, Double_t psi3);
         // analysis details
         Bool_t                  CorrectRho(Double_t psi2, Double_t psi3);
         // event and track selection, also used by AliAnalyisTaskJetFlow
@@ -163,6 +168,7 @@ class AliAnalysisTaskRhoVnModulation : public AliAnalysisTaskEmcalJet
         TH1F*                   GetResolutionFromOuptutFile(detectorType detector, Int_t h = 2, TArrayD* c = 0x0);
         TH1F*                   CorrectForResolutionDiff(TH1F* v, detectorType detector, TArrayD* cen, Int_t c, Int_t h = 2);
         TH1F*                   CorrectForResolutionInt(TH1F* v, detectorType detector, TArrayD* cen, Int_t h = 2);
+        TH1F*                   GetDifferentialQC(TProfile* refCumulants, TProfile* diffCumlants, TArrayD* ptBins, Int_t h);
     private:
         // analysis flags and settings
         Int_t                   fDebug;                 // debug level (0 none, 1 fcn calls, 2 verbose)
@@ -170,6 +176,7 @@ class AliAnalysisTaskRhoVnModulation : public AliAnalysisTaskEmcalJet
         Bool_t                  fFillQAHistograms;      // fill qa histograms
         Int_t                   fReduceBinsXByFactor;   // reduce the bins on x-axis of histo's by this integer
         Int_t                   fReduceBinsYByFactor;   // reduce the bins on y-axis of histo's by this integer
+        Bool_t                  fNoEventWeightsForQC;   // don't store event weights for qc analysis
         TArrayI*                fCentralityClasses;     //-> centrality classes (maximum 10)
         TArrayD*                fPtBinsHybrids;         //-> pt bins for hybrid track vn anaysis
         TArrayD*                fPtBinsJets;            //-> pt bins for jet vn analysis
@@ -181,6 +188,7 @@ class AliAnalysisTaskRhoVnModulation : public AliAnalysisTaskEmcalJet
         Int_t                   fNAcceptedTracks;       //! number of accepted tracks
         Int_t                   fNAcceptedTracksQCn;    //! accepted tracks for QCn
         fitModulationType       fFitModulationType;     // fit modulation type
+        qcRecovery              fQCRecovery;            // recovery type for e-by-e qc method
         Bool_t                  fUsePtWeight;           // use dptdphi instead of dndphi
         detectorType            fDetectorType;          // type of detector used for modulation fit
         TString                 fFitModulationOptions;  // fit options for modulation fit
@@ -236,6 +244,7 @@ class AliAnalysisTaskRhoVnModulation : public AliAnalysisTaskEmcalJet
         TProfile*               fProfV3Resolution[10];  //! resolution parameters for v3
         // qa histograms for accepted pico tracks
         TH1F*                   fHistPicoTrackPt[10];    //! pt of all charged tracks
+        TH1F*                   fHistPicoTrackMult[10];  //! multiplicity of accepted pico tracks
         TH2F*                   fHistPicoCat1[10];       //! pico tracks spd hit and refit
         TH2F*                   fHistPicoCat2[10];       //! pico tracks wo spd hit w refit, constrained
         TH2F*                   fHistPicoCat3[10];       //! pico tracks wo spd hit wo refit, constrained
@@ -306,7 +315,7 @@ class AliAnalysisTaskRhoVnModulation : public AliAnalysisTaskEmcalJet
         AliAnalysisTaskRhoVnModulation(const AliAnalysisTaskRhoVnModulation&);                  // not implemented
         AliAnalysisTaskRhoVnModulation& operator=(const AliAnalysisTaskRhoVnModulation&);       // not implemented
 
-        ClassDef(AliAnalysisTaskRhoVnModulation, 13);
+        ClassDef(AliAnalysisTaskRhoVnModulation, 14);
 };
 
 #endif