typo fix in addtask and minor update for dphi dpt fits
authorrbertens <rbertens@cern.ch>
Wed, 23 Jul 2014 10:28:52 +0000 (12:28 +0200)
committerrbertens <rbertens@cern.ch>
Wed, 23 Jul 2014 10:29:22 +0000 (12:29 +0200)
PWG/FLOW/Tasks/AliJetFlowTools.cxx
PWG/FLOW/Tasks/AliJetFlowTools.h
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetV2.cxx
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetV2.h
PWGJE/EMCALJetTasks/macros/AddTaskJetV2.C

index 4474ad1..c71b5eb 100644 (file)
@@ -78,6 +78,7 @@ AliJetFlowTools::AliJetFlowTools() :
     fResponseMaker      (new AliAnaChargedJetResponseMaker()),
     fRMS                (kTRUE),
     fSymmRMS            (kTRUE),
+    fRho0               (kFALSE),
     fPower              (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
     fSaveFull           (kTRUE),
     fActiveString       (""),
@@ -882,6 +883,7 @@ Bool_t AliJetFlowTools::PrepareForUnfolding()
     }
     // 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;
@@ -949,7 +951,11 @@ Bool_t AliJetFlowTools::PrepareForUnfolding()
     }
     // 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());
@@ -1031,8 +1037,9 @@ Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) {
         // 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());
@@ -1049,7 +1056,11 @@ Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) {
     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());
@@ -1542,6 +1553,10 @@ void AliJetFlowTools::Style(TH1* h, EColor col, histoType type, Bool_t legacy)
             h->SetMarkerStyle(8);
             h->SetMarkerSize(1);
        } break;
+       case kDeltaPhi : {
+            h->GetYaxis()->SetTitle("[counts]");
+            h->GetXaxis()->SetTitle("#Delta #phi");
+       }
        default : break;
     }
 }
@@ -3819,10 +3834,16 @@ void AliJetFlowTools::MakeAU() {
     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;
@@ -3870,7 +3891,8 @@ void AliJetFlowTools::MakeAU() {
             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]");
@@ -3900,9 +3922,11 @@ void AliJetFlowTools::MakeAU() {
         
         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);
         }
@@ -3912,13 +3936,28 @@ void AliJetFlowTools::MakeAU() {
     }
     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) {
index 26aa04d..39b190a 100644 (file)
@@ -60,6 +60,7 @@ class AliJetFlowTools {
             kBar,                       // default style for bar histogram
             kRatio,                     // default style for ratio
             kV2,                        // default style for v2
+            kDeltaPhi,                  // default for delta phi
             kEmpty };                   // default style
         // setters, interface to the class
         void            SetSaveFull(Bool_t b)           {fSaveFull              = b;}
@@ -138,6 +139,7 @@ class AliJetFlowTools {
         void            SetWeightFunction(TF1* w)               {fResponseMaker->SetRMMergeWeightFunction(w);}
         void            SetRMS(Bool_t r)                        {fRMS                   = r;}
         void            SetSymmRMS(Bool_t r)                    {fSymmRMS               = r; fRMS               = r;}
+        void            SetRho0(Bool_t r)                       {fRho0                  = r;}
         void            Make();
         void            MakeAU();       // test function, use with caution (09012014)
         void            Finish() {
@@ -389,6 +391,7 @@ TLatex* tex = new TLatex(xmin, ymax, string.Data());
         AliAnaChargedJetResponseMaker*  fResponseMaker; // utility object
         Bool_t                  fRMS;                   // systematic method
         Bool_t                  fSymmRMS;               // symmetric systematic method
+        Bool_t                  fRho0;                  // use the result obtained with the 'classic' fixed rho
         TF1*                    fPower;                 // smoothening fit
         Bool_t                  fSaveFull;              // save all generated histograms to file
         TString                 fActiveString;          // identifier of active output
index 0175f05..625e72a 100644 (file)
@@ -707,6 +707,7 @@ Bool_t AliAnalysisTaskJetV2::Run()
 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);
@@ -770,6 +771,7 @@ void AliAnalysisTaskJetV2::CalculateEventPlaneTPC(Double_t* tpc)
 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);
@@ -1181,7 +1183,7 @@ Bool_t AliAnalysisTaskJetV2::CorrectRho(Double_t psi2, Double_t psi3)
     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) {
index ec952fc..59de265 100644 (file)
@@ -30,11 +30,11 @@ class AliAnalysisTaskJetV2 : public AliAnalysisTaskEmcalJet {
          // 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();
@@ -77,6 +77,7 @@ class AliAnalysisTaskJetV2 : public AliAnalysisTaskEmcalJet {
         }
         /* 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");
index 5c5bae9..ca53819 100644 (file)
@@ -117,7 +117,7 @@ AliAnalysisTaskJetV2* AddTaskJetV2(
   (LHC10h) ? jetTask->SetExpectedSemiGoodRuns(0x0) : jetTask->SetExpectedSemiGoodRuns(new TArrayI(sizeof(semiGoodRuns)/sizeof(semiGoodRuns[0]), semiGoodRuns));
 
   // and if 10h, pass this info to the task so acceptance isn't changed
-  if(LHC10h) task->SetCollisionType(AliAnalysisTaskJetV2::kPbPb10h);
+  if(LHC10h) jetTask->SetCollisionType(AliAnalysisTaskJetV2::kPbPb10h);
 
   //-------------------------------------------------------
   // Final settings, pass to manager and set the containers