flow package documentation cleanup, dphi-dpt unfolding, bugfix: fixed crash when...
authorrbertens <rbertens@cern.ch>
Mon, 13 Jan 2014 15:54:33 +0000 (16:54 +0100)
committerrbertens <rbertens@cern.ch>
Mon, 13 Jan 2014 15:54:33 +0000 (16:54 +0100)
20 files changed:
PWG/FLOW/Tasks/AliJetFlowTools.cxx
PWG/FLOW/Tasks/AliJetFlowTools.h
PWGCF/FLOW/Documentation/FLOW.tex
PWGCF/FLOW/Documentation/chapters/FitQDistribution.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/FlowEvent.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/FrontMatter.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/GeneratingFunctionCumulants.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/Introduction.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/LeeYangZeros.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/LeeYangZerosEP.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/Methods.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/MonteCarlo.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/Program.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/QVectorCumulants.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/QuickStart.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/ScalarProduct.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/Summary.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/appendix1.tex [deleted file]
PWGCF/FLOW/Documentation/chapters/bibliography.tex [deleted file]
PWGCF/FLOW/Documentation/figs/flowChart.png [new file with mode: 0644]

index e2ed0e3..7a09f57 100644 (file)
@@ -112,7 +112,9 @@ AliJetFlowTools::AliJetFlowTools() :
     fUseDetectorResponse(kTRUE),
     fUseDptResponse     (kTRUE),
     fTrainPower         (kTRUE),
-    fInOutUnfolding     (kTRUE),
+    fDphiUnfolding      (kTRUE),
+    fDphiDptUnfolding   (kFALSE),
+    fExLJDpt            (kTRUE),
     fRMSSpectrumIn      (0x0),
     fRMSSpectrumOut     (0x0),
     fRMSRatio           (0x0),
@@ -127,12 +129,18 @@ AliJetFlowTools::AliJetFlowTools() :
     fDptOut             (0x0),
     fFullResponseIn     (0x0),
     fFullResponseOut    (0x0) { // class constructor
-    // create response maker weight function (tuned to PYTHIA spectrum)
-    fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0 ,200));
+        // create response maker weight function (tuned to PYTHIA spectrum)
+        fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0, 200));
+        for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
 }
 //_____________________________________________________________________________
 void AliJetFlowTools::Make() {
     // core function of the class
+    if(fDphiDptUnfolding) {
+        // to extract the yield as function of Dphi, Dpt - experimental
+        MakeAU();
+        return;
+    }
     // 1) rebin the raw output of the jet task to the desired binnings
     // 2) calls the unfolding routine
     // 3) writes output to file
@@ -206,58 +214,14 @@ void AliJetFlowTools::Make() {
     fActiveDir->cd();                   // select active dir
     TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
     dirIn->cd();                        // select inplane subdir
-    // select the unfolding method
-    switch (fUnfoldingAlgorithm) {
-        case kChi2 : {
-            unfoldedJetSpectrumIn = UnfoldSpectrumChi2(       // do the inplane unfolding
-                measuredJetSpectrumIn,
-                resizedResponseIn,
-                kinematicEfficiencyIn,
-                measuredJetSpectrumTrueBinsIn,
-                TString("in"),
-                jetFindingEfficiency);
-            printf(" > Spectrum (in plane) unfolded using kChi2 unfolding < \n");
-        } break;
-        case kBayesian : {
-            unfoldedJetSpectrumIn = UnfoldSpectrumBayesian(       // do the inplane unfolding
-                measuredJetSpectrumIn,
-                resizedResponseIn,
-                kinematicEfficiencyIn,
-                measuredJetSpectrumTrueBinsIn,
-                TString("in"),
-                jetFindingEfficiency);
-            printf(" > Spectrum (in plane) unfolded using kBayesian unfolding < \n");
-        } break;
-        case kBayesianAli : {
-            unfoldedJetSpectrumIn = UnfoldSpectrumBayesianAli(       // do the inplane unfolding
-                measuredJetSpectrumIn,
-                resizedResponseIn,
-                kinematicEfficiencyIn,
-                measuredJetSpectrumTrueBinsIn,
-                TString("in"),
-                jetFindingEfficiency);
-            printf(" > Spectrum (in plane) unfolded using kBayesianAli unfolding < \n");
-        } break;
-        case kSVD : {
-            unfoldedJetSpectrumIn = UnfoldSpectrumSVD(       // do the inplane unfolding
-                measuredJetSpectrumIn,
-                resizedResponseIn,
-                kinematicEfficiencyIn,
-                measuredJetSpectrumTrueBinsIn,
-                TString("in"),
-                jetFindingEfficiency);
-            printf(" > Spectrum (in plane) unfolded using kSVD unfolding < \n");
-        } break;
-        case kNone : {  // do nothing 
-            resizedResponseIn->SetNameTitle("measuredSpectrumIn", "measured spectrum, in plane");
-            unfoldedJetSpectrumIn = ProtectHeap(measuredJetSpectrumIn, kTRUE, TString("in"));
-        } break;
-        
-        default : {
-            printf(" > Selected unfolding method is not implemented yet ! \n");
-            return;
-        }
-    }
+    // do the inplane unfolding
+    unfoldedJetSpectrumIn = UnfoldWrapper(
+        measuredJetSpectrumIn,
+        resizedResponseIn,
+        kinematicEfficiencyIn,
+        measuredJetSpectrumTrueBinsIn,
+        TString("in"),
+        jetFindingEfficiency);
     resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
     resizedResponseIn->SetXTitle("p_{T}^{true} [GeV/c]");
     resizedResponseIn->SetYTitle("p_{T}^{rec} [GeV/c]");
@@ -281,59 +245,17 @@ void AliJetFlowTools::Make() {
         fFullResponseIn->Write();
     }
     fActiveDir->cd();
-    if(fInOutUnfolding) {
+    if(fDphiUnfolding) {
         TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
         dirOut->cd();
-        switch (fUnfoldingAlgorithm) {
-            case kChi2 : {
-                unfoldedJetSpectrumOut = UnfoldSpectrumChi2(
-                    measuredJetSpectrumOut,
-                    resizedResponseOut,
-                    kinematicEfficiencyOut,
-                    measuredJetSpectrumTrueBinsOut,
-                    TString("out"),
-                    jetFindingEfficiency);
-                printf(" > Spectrum (out of plane) unfolded using kChi2 < \n");
-            } break;
-            case kBayesian : {
-                unfoldedJetSpectrumOut = UnfoldSpectrumBayesian(
-                    measuredJetSpectrumOut,
-                    resizedResponseOut,
-                    kinematicEfficiencyOut,
-                    measuredJetSpectrumTrueBinsOut,
-                    TString("out"),
-                    jetFindingEfficiency);
-                printf(" > Spectrum (out of plane) unfolded using kBayesian < \n");
-            } break;
-            case kBayesianAli : {
-                unfoldedJetSpectrumOut = UnfoldSpectrumBayesianAli(
+        // do the out of plane unfolding
+        unfoldedJetSpectrumOut = UnfoldWrapper(
                     measuredJetSpectrumOut,
                     resizedResponseOut,
                     kinematicEfficiencyOut,
                     measuredJetSpectrumTrueBinsOut,
                     TString("out"),
                     jetFindingEfficiency);
-                printf(" > Spectrum (out of plane) unfolded using kBayesianAli < \n");
-            } break;
-            case kSVD : {
-                unfoldedJetSpectrumOut = UnfoldSpectrumSVD(
-                    measuredJetSpectrumOut,
-                    resizedResponseOut,
-                    kinematicEfficiencyOut,
-                    measuredJetSpectrumTrueBinsOut,
-                    TString("out"),
-                    jetFindingEfficiency);
-                printf(" > Spectrum (out of plane) unfolded using kSVD < \n");
-            } break;
-            case kNone : {  // do nothing
-                resizedResponseOut->SetNameTitle("measuredSpectrumOut", "measured spectrum, out plane");
-                unfoldedJetSpectrumOut = ProtectHeap(measuredJetSpectrumOut, kTRUE, TString("out"));
-            } break;
-            default : {
-                printf(" > Selected unfolding method is not implemented yet ! \n");
-                return;
-            }
-        }
         resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
         resizedResponseOut->SetXTitle("p_{T}^{true} [GeV/c]");
         resizedResponseOut->SetYTitle("p_{T}^{rec} [GeV/c]");
@@ -402,16 +324,57 @@ void AliJetFlowTools::Make() {
                 v2->Write();
             }
         }
-    }   // end of if(fInOutUnfolding)
+    }   // end of if(fDphiUnfolding)
     fDeltaPtDeltaPhi->Write();
     fJetPtDeltaPhi->Write();
     // save the current state of the unfolding object
     SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, unfoldedJetSpectrumOut ? kTRUE : kFALSE);
 }
 //_____________________________________________________________________________
+TH1D* AliJetFlowTools::UnfoldWrapper(
+        const TH1D* measuredJetSpectrum,        // truncated raw jets (same binning as pt rec of response) 
+        const TH2D* resizedResponse,            // response matrix
+        const TH1D* kinematicEfficiency,        // kinematic efficiency
+        const TH1D* measuredJetSpectrumTrueBins,        // unfolding template: same binning is pt gen of response
+        const TString suffix,                   // suffix (in or out of plane)
+        const TH1D* jetFindingEfficiency)       // jet finding efficiency
+{
+    // wrapper function to call specific unfolding routine
+    TH1D* (AliJetFlowTools::*myFunction)(const TH1D*, const TH2D*, const TH1D*, const TH1D*, const TString, const TH1D*);
+    // initialize functon pointer
+    switch (fUnfoldingAlgorithm) {
+        case kChi2 : {
+            myFunction = &AliJetFlowTools::UnfoldSpectrumChi2;
+            break; 
+        }
+        case kBayesian : {
+            myFunction = &AliJetFlowTools::UnfoldSpectrumBayesian;
+            break;
+        }
+        case kBayesianAli : {
+            myFunction = &AliJetFlowTools::UnfoldSpectrumBayesianAli;
+            break;
+        }
+        case kSVD : {
+            myFunction = &AliJetFlowTools::UnfoldSpectrumSVD;
+            break;
+        }
+        case kNone : {
+            TH1D* clone((TH1D*)resizedResponse->Clone("clone"));
+            clone->SetNameTitle(Form("measuredSpectrum_%s", suffix.Data()), Form("measured spectrum %s", suffix.Data()));
+            return clone;
+        }
+        default : {
+            return 0x0;
+        }
+    }
+    // do the actual unfolding with the selected function
+    return (this->*myFunction)( measuredJetSpectrum, resizedResponse, kinematicEfficiency, measuredJetSpectrumTrueBins, suffix, jetFindingEfficiency);
+} 
+//_____________________________________________________________________________
 TH1D* AliJetFlowTools::UnfoldSpectrumChi2(
         const TH1D* measuredJetSpectrum,      // truncated raw jets (same binning as pt rec of response) 
-        const TH2D* resizedResponse,           // response matrix
+        const TH2D* resizedResponse,          // response matrix
         const TH1D* kinematicEfficiency,      // kinematic efficiency
         const TH1D* measuredJetSpectrumTrueBins,        // unfolding template: same binning is pt gen of response
         const TString suffix,                 // suffix (in or out of plane)
@@ -891,7 +854,7 @@ Bool_t AliJetFlowTools::PrepareForUnfolding()
         printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
         return kFALSE;
     }
-    if(!fRMSSpectrumIn && fInOutUnfolding) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
+    if(!fRMSSpectrumIn && fDphiUnfolding) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
                           // procedures, these profiles will be nonsensical, user is responsible
         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
@@ -910,7 +873,7 @@ Bool_t AliJetFlowTools::PrepareForUnfolding()
     }
     fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
     // in plane spectrum
-    if(!fInOutUnfolding) {
+    if(!fDphiUnfolding) {
         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40, "e");
         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40, "e");
     } else {
@@ -952,7 +915,8 @@ Bool_t AliJetFlowTools::PrepareForUnfolding()
         if(normalizeToFullSpectrum) fEventCount = -1;
     }
     // extract the delta pt matrices
-    TString deltaptName(Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityBin));
+    TString deltaptName("");
+    deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityBin) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityBin);
     fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
     if(!fDeltaPtDeltaPhi) {
         printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
@@ -962,7 +926,7 @@ Bool_t AliJetFlowTools::PrepareForUnfolding()
     }
     fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
     // in plane delta pt distribution
-    if(!fInOutUnfolding) {
+    if(!fDphiUnfolding) {
         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40, "e");
         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40, "e");
     } else {
@@ -1007,6 +971,64 @@ Bool_t AliJetFlowTools::PrepareForUnfolding()
     return kTRUE;
 }
 //_____________________________________________________________________________
+Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) {
+    // prepare for unfoldingUA - more robust method to extract input spectra from file
+    // will replace PrepareForUnfolding eventually (09012014)
+    if(!fInputList) {
+        printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
+        return kFALSE;
+    }
+    if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
+    // check if the pt bin for true and rec have been set
+    if(!fBinsTrue || !fBinsRec) {
+        printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
+        return kFALSE;
+    }
+    if(!fTrainPower) {
+        // 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
+    TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityBin));
+    fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
+    if(!fJetPtDeltaPhi) {
+        printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
+        return kFALSE;
+    }
+    fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
+    // in plane spectrum
+    fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), low, up, "e");
+    // extract the delta pt matrices
+    TString deltaptName("");
+    deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityBin) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityBin);
+    fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
+    if(!fDeltaPtDeltaPhi) {
+        printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
+        printf(" > may be ok, depending no what you want to do < \n");
+        fRefreshInput = kTRUE;
+        return kTRUE;
+    }
+    fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
+    // in plane delta pt distribution
+    fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), low, up, "e");
+    // create a rec - true smeared response matrix
+    TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
+    for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
+        Bool_t skip = kFALSE;
+        for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
+            (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
+            if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
+        }
+    }
+    fDptIn = new TH2D(*rfIn);
+    fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
+    fDptIn->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
+    fDptIn->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
+    fDptIn = ProtectHeap(fDptIn);
+    
+    return kTRUE;
+}
+//_____________________________________________________________________________
 TH1D* AliJetFlowTools::GetPrior(
         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins 
         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
@@ -1234,12 +1256,11 @@ TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min
     if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
     TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
     temp->Sumw2();      // if already called on the original, this will give off a warning but do nothing
-    TFitResultPtr r = temp->Fit(function, "WLIS", "", min, max);
+    TFitResultPtr r = temp->Fit(function, "VIS", "", min, max);
     if((int)r == 0) {   // MINUIT status
         for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
             if(temp->GetBinCenter(i) > start) {     // from this pt value use extrapolation
-                if(counts) temp->SetBinContent(i, (int)function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
-                else temp->SetBinContent(i, function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
+                (counts) ? temp->SetBinContent(i, (int)(function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i))) : temp->SetBinContent(i, function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
                 if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
             }
         }
@@ -1384,7 +1405,7 @@ void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type)
     }
 }
 //_____________________________________________________________________________
-void AliJetFlowTools::PostProcess(TString def, TString in, TString out, Int_t columns) const
+void AliJetFlowTools::PostProcess(TString def, Int_t columns, TString in, TString out) const
 {
    // go through the output file and perform post processing routines
    // can either be performed in one go with the unfolding, or at a later stage
@@ -1404,22 +1425,22 @@ void AliJetFlowTools::PostProcess(TString def, TString in, TString out, Int_t co
    // prepare necessary canvasses
    TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
    TCanvas* canvasOut(0x0);
-   if(fInOutUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
+   if(fDphiUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
    TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
-   if(fInOutUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
+   if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
    TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn")); 
    TCanvas* canvasSpectraOut(0x0);
-   if(fInOutUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
+   if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
    TCanvas* canvasRatio(0x0);
-   if(fInOutUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
+   if(fDphiUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
    TCanvas* canvasV2(0x0);
-   if(fInOutUnfolding) canvasV2 = new TCanvas("V2", "V2");
+   if(fDphiUnfolding) canvasV2 = new TCanvas("V2", "V2");
    TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
    TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
    TCanvas* canvasMasterOut(0x0);
-   if(fInOutUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
-   (fInOutUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
+   if(fDphiUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
+   (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
    TDirectoryFile* defDir(0x0);
    
    // get an estimate of the number of outputs and find the default set
@@ -1466,6 +1487,102 @@ void AliJetFlowTools::PostProcess(TString def, TString in, TString out, Int_t co
        tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
        tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
        j++;
+       if(!tempIn) {    // attempt to get MakeAU output
+           TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
+           TCanvas* tempPearson(new TCanvas(Form("pearson_%s", dirName.Data()), Form("pearson_%s", dirName.Data())));
+           tempPearson->Divide(4,2);
+           TCanvas* tempRatio(new TCanvas(Form("ratio_%s", dirName.Data()), Form("ratio_%s", dirName.Data())));
+           tempRatio->Divide(4,2);
+           TCanvas* tempResult(new TCanvas(Form("result_%s", dirName.Data()), Form("result_%s", dirName.Data())));
+           tempResult->Divide(4,2);
+           for(Int_t q(0); q < 8; q++) {
+               tempIn = (TDirectoryFile*)tempDir->Get(Form("%s___%s", stringArray[q].Data(), dirName.Data()));
+               if(tempIn) {
+                    // to see if the unfolding converged try to extract the pearson coefficients
+                    TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
+                    if(pIn) {
+                       printf(" - %s in plane converged \n", dirName.Data());
+                        tempPearson->cd(1+q);
+                        Style(gPad, "PEARSON");
+                        pIn->DrawCopy("colz");
+                        TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
+                        if(rIn) {
+                            Style(rIn);
+                            printf(" > found RatioRefoldedMeasured < \n");
+                            tempRatio->cd(q+1);
+                            rIn->SetFillColor(kRed);
+                            rIn->Draw("ap");
+                        }
+                        TH1D* dvector((TH1D*)tempIn->Get("dVector"));
+                        TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
+                        TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
+                        TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
+                        if(dvector && avalue && rm && eff) {
+                            Style(dvector);
+                            Style(avalue);
+                            Style(rm);
+                            Style(eff);
+                            canvasMISC->cd(1);
+                            Style(gPad, "SPECTRUM");
+                            dvector->DrawCopy();
+                            canvasMISC->cd(2);
+                            Style(gPad, "SPECTRUM");
+                            avalue->DrawCopy();
+                            canvasMISC->cd(3);
+                            Style(gPad, "PEARSON");
+                            rm->DrawCopy("colz");
+                            canvasMISC->cd(4);
+                            eff->DrawCopy();
+                        } else if(rm && eff) {
+                            Style(rm);
+                            Style(eff);
+                            canvasMISC->cd(3);
+                            Style(gPad, "PEARSON");
+                            rm->DrawCopy("colz");
+                            canvasMISC->cd(4);
+                            eff->DrawCopy();
+                        }
+                    }
+                    TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
+                    TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
+                    TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
+                    if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
+                        if(defaultUnfoldedJetSpectrumIn) {
+                            Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
+                            TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
+                            temp->Divide(unfoldedSpectrum);
+                            temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
+                            temp->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+                            temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
+                            canvasMasterIn->cd(j);
+                            temp->GetYaxis()->SetRangeUser(0., 2);
+                            temp->DrawCopy();
+                        }
+                        TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
+                        tempResult->cd(q+1);
+                        Style(gPad);
+                        Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
+                        unfoldedSpectrum->DrawCopy();
+                        Style(inputSpectrum, kGreen, kMeasuredSpectrum);
+                        inputSpectrum->DrawCopy("same");
+                        Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
+                        refoldedSpectrum->DrawCopy("same");
+                        TLegend* l(AddLegend(gPad));
+                        Style(l);
+                        if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
+                            Float_t chi(fitStatus->GetBinContent(1));
+                            Float_t pen(fitStatus->GetBinContent(2));
+                            Int_t dof((int)fitStatus->GetBinContent(3));
+                            Float_t beta(fitStatus->GetBinContent(4));
+                            l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
+                        } else if (fitStatus) { // only available in SVD fit
+                            Int_t reg((int)fitStatus->GetBinContent(1));
+                            l->AddEntry((TObject*)0, Form("REG %i", reg), "");
+                        }
+                    }
+               }
+           }
+       }
        if(tempIn) { 
            // to see if the unfolding converged try to extract the pearson coefficients
            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
@@ -1758,7 +1875,11 @@ TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t a
     gr->GetXaxis()->SetTitle("p_{T} [GeV/c]");
     Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
     TH1* dud((TH1*)h1->Clone("dud"));
-    if(!dud->Divide(h1, h2, 1., 1., "B")) {
+    dud->Sumw2();
+    h1->Sumw2();
+    h2->Sumw2();
+    if(!dud->Divide(h1, h2)) {
+        printf(" > ROOT failed to divide two histogams, dividing manually \n < ");
         for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
             binCent = h1->GetXaxis()->GetBinCenter(i);
             if(xmax > 0. && binCent > xmax) continue;
@@ -1766,7 +1887,7 @@ TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t a
             binWidth = h1->GetXaxis()->GetBinWidth(i);
             if(h2->GetBinContent(j) > 0.) {
                 ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
-            /* original propagation of uncertainty
+            /* original propagation of uncertainty changed 08012014
             Double_t A = 1./h2->GetBinContent(j)*h1->GetBinError(i);
             Double_t B = 0.;
             if(h2->GetBinError(j)>0.) {
@@ -1777,11 +1898,12 @@ TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t a
                 Double_t B = h2->GetBinError(i)/h2->GetBinContent(i);
                 error2 = ratio*ratio*A*A+ratio*ratio*B*B;
                 if(error2 > 0 ) error2 = TMath::Sqrt(error2);
-                gr->SetPoint(gr->GetN(),binCent,ratio);
-                gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
+                gr->SetPoint(gr->GetN(), binCent, ratio);
+                gr->SetPointError(gr->GetN()-1, 0.5*binWidth, error2);
             }
         }
     } else {
+        printf( " > using ROOT to divide two histograms < \n");
         for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
             binCent = dud->GetXaxis()->GetBinCenter(i);
             if(xmax > 0. && binCent > xmax) continue;
@@ -2026,3 +2148,121 @@ TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, T
     return p;
 }
 //_____________________________________________________________________________
+void AliJetFlowTools::MakeAU() {
+    // === azimuthal unfolding ===
+    // 
+    // unfolds the spectrum in delta phi bins, extracts the yield per bin, and does a fit
+    // in transverse momentum and azimuthal correlation space to extract v2 params
+    // settings are equal to the ones used for 'Make()'
+    //
+    // basic steps that are followed:
+    // 1) rebin the raw output of the jet task to the desired binnings
+    // 2) calls the unfolding routine
+    // 3) writes output to file
+    // can be repeated multiple times with different configurations
+
+    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[6];
+    for(Int_t i(0); i < 6; i++) dPtdPhi[i] = new TH1D(Form("dPtdPhi_%i", i), Form("dPtdPhi_%i", i), 8, 0, TMath::Pi());
+
+    for(Int_t i(0); i < 8; i++) {
+        // 1) manipulation of input histograms
+        // check if the input variables are present
+        if(!PrepareForUnfolding(low[i], up[i])) return;
+        // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
+        //     parts of the spectrum can end up in over or underflow bins
+        TH1D* measuredJetSpectrumIn  = RebinTH1D(fSpectrumIn, fBinsRec, Form("resized_%s", stringArray[i].Data()), kFALSE);
+
+        // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
+        // the template will be used as a prior for the chi2 unfolding
+        TH1D* measuredJetSpectrumTrueBinsIn  = RebinTH1D(fSpectrumIn, fBinsTrue, stringArray[i], kFALSE);
+
+        // get the full response matrix from the dpt and the detector response
+        fDetectorResponse = NormalizeTH2D(fDetectorResponse);
+        // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
+        // so that unfolding should return the initial spectrum
+        if(fUseDptResponse && fUseDetectorResponse) fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
+        else if (fUseDptResponse && !fUseDetectorResponse) fFullResponseIn = fDptIn;
+        else if (!fUseDptResponse && fUseDetectorResponse) fFullResponseIn = fDetectorResponse;
+        else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) return;
+        // normalize each slide of the response to one
+        NormalizeTH2D(fFullResponseIn);
+        // resize to desired binning scheme
+        TH2D* resizedResponseIn  = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, stringArray[i]);
+        // get the kinematic efficiency
+        TH1D* kinematicEfficiencyIn  = resizedResponseIn->ProjectionX();
+        kinematicEfficiencyIn->SetNameTitle(Form("kin_eff_%s", stringArray[i].Data()), Form("kin_eff_%s", stringArray[i].Data()));
+        // suppress the errors 
+        for(Int_t j(0); j < kinematicEfficiencyIn->GetXaxis()->GetNbins(); j++) kinematicEfficiencyIn->SetBinError(1+j, 0.);
+        TH1D* jetFindingEfficiency(0x0);
+        if(fJetFindingEff) {
+            jetFindingEfficiency = ProtectHeap(fJetFindingEff);
+            jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
+            jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
+        }
+        // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
+        TH1D* unfoldedJetSpectrumIn(0x0);
+        fActiveDir->cd();                   // select active dir
+        TDirectoryFile* dirIn = new TDirectoryFile(Form("%s___%s", stringArray[i].Data(), fActiveString.Data()), Form("%s___%s", stringArray[i].Data(), fActiveString.Data()));
+        dirIn->cd();                        // select inplane subdir
+        // select the unfolding method
+        unfoldedJetSpectrumIn = UnfoldWrapper(
+            measuredJetSpectrumIn,
+            resizedResponseIn,
+            kinematicEfficiencyIn,
+            measuredJetSpectrumTrueBinsIn,
+            TString("dPtdPhiUnfolding"),
+            jetFindingEfficiency);
+        if(i==5) {
+            resizedResponseIn->SetNameTitle(Form("ResponseMatrix_%s", stringArray[i].Data()), Form("response matrix %s", stringArray[i].Data()));
+            resizedResponseIn->SetXTitle("p_{T}^{true} [GeV/c]");
+            resizedResponseIn->SetYTitle("p_{T}^{rec} [GeV/c]");
+            resizedResponseIn = ProtectHeap(resizedResponseIn);
+            resizedResponseIn->Write();
+            kinematicEfficiencyIn->SetNameTitle(Form("KinematicEfficiency_%s", stringArray[i].Data()), Form("Kinematic efficiency, %s", stringArray[i].Data()));
+            kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
+            kinematicEfficiencyIn->Write();
+            fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
+            fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
+            fDetectorResponse->Write();
+            // optional histograms
+            if(fSaveFull) {
+                fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", Form("[INPUT] Jet spectrum, %s", stringArray[i].Data()));
+                fSpectrumIn->Write();
+                fDptInDist->SetNameTitle("[ORIG]DeltaPt", Form("#delta p_{T} distribution, %s", stringArray[i].Data()));
+                fDptInDist->Write();
+                fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix", Form("#delta p_{T} matrix, %s", stringArray[i].Data()));
+                fDptIn->Write();
+                fFullResponseIn->SetNameTitle("ResponseMatrix", Form("Response matrix, %s", stringArray[i].Data()));
+                fFullResponseIn->Write();
+            }
+        }
+        fActiveDir->cd();
+        fDeltaPtDeltaPhi->Write();
+        fJetPtDeltaPhi->Write();
+        
+        TH1D* dud(ProtectHeap(unfoldedJetSpectrumIn, kTRUE, stringArray[i]));;
+        Double_t integralError(0);
+        for(Int_t j(0); j < 6; j++) {
+            // get the integrated 
+           Double_t integral(dud->IntegralAndError(2*j+1, 2*j+3, integralError));
+           dPtdPhi[j]->SetBinContent(i+1, integral);
+           dPtdPhi[j]->SetBinError(i+1, integralError);
+        }
+        dud->Write();
+        // save the current state of the unfolding object
+        SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, kFALSE);
+    }
+    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++) {
+        dPtdPhi[i]->Fit(fourier, "VI");
+        v2->SetBinContent(1+i, fourier->GetParameter(1));
+        v2->SetBinError(1+i, fourier->GetParError(1));
+        dPtdPhi[i]->Write();
+    }
+    v2->Write();
+}
+//_____________________________________________________________________________
index f220cef..66a0a65 100644 (file)
@@ -109,8 +109,12 @@ class AliJetFlowTools {
         void            SetUseDetectorResponse(Bool_t r)        {fUseDetectorResponse   = r;}
         void            SetUseDptResponse(Bool_t r)             {fUseDptResponse        = r;}
         void            SetTrainPowerFit(Bool_t t)              {fTrainPower            = t;}
-        void            SetUnfoldInOutPlane(Bool_t i)           {fInOutUnfolding        = i;}
+        void            SetDphiUnfolding(Bool_t i)              {fDphiUnfolding         = i;}
+        void            SetDphiDptUnfolding(Bool_t i)           {fDphiDptUnfolding      = i;}
+        void            SetExLJDpt(Bool_t i)                    {fExLJDpt               = i;}
+        void            SetWeightFunction(TF1* w)               {fResponseMaker->SetRMMergeWeightFunction(w);}
         void            Make();
+        void            MakeAU();       // test function, use with caution (09012014)
         void            Finish() {
             fOutputFile->cd();
             if(fRMSSpectrumIn)  fRMSSpectrumIn->Write();
@@ -119,9 +123,9 @@ class AliJetFlowTools {
             fOutputFile->Close();}
         void            PostProcess(
                 TString def,
+                Int_t columns = 4,
                 TString in = "UnfoldedSpectra.root", 
-                TString out = "ProcessedSpectra.root", 
-                Int_t columns = 4) const;
+                TString out = "ProcessedSpectra.root") const;
         Bool_t          SetRawInput (
                 TH2D* detectorResponse, // detector response matrix
                 TH1D* jetPtIn,          // in plane jet spectrum
@@ -161,12 +165,19 @@ class AliJetFlowTools {
         static void     SetDebug(Int_t d)               {AliUnfolding::SetDebug(d);}
     private:
         Bool_t          PrepareForUnfolding(); 
+        Bool_t          PrepareForUnfolding(Int_t low, Int_t up);
         TH1D*           GetPrior(                       const TH1D* measuredJetSpectrum,
                                                         const TH2D* resizedResponse,
                                                         const TH1D* kinematicEfficiency,
                                                         const TH1D* measuredJetSpectrumTrueBins,
                                                         const TString suffix,
                                                         const TH1D* jetFindingEfficiency);
+        TH1D*           UnfoldWrapper(                  const TH1D* measuredJetSpectrum, 
+                                                        const TH2D* resizedResponse,
+                                                        const TH1D* kinematicEfficiency,
+                                                        const TH1D* measuredJetSpectrumTrueBins,
+                                                        const TString suffix,
+                                                        const TH1D* jetFindingEfficiency = 0x0);
         TH1D*           UnfoldSpectrumChi2(             const TH1D* measuredJetSpectrum, 
                                                         const TH2D* resizedResponse,
                                                         const TH1D* kinematicEfficiency,
@@ -236,15 +247,16 @@ class AliJetFlowTools {
         Float_t                 fFitStart;              // from this value, use smoothening
         Bool_t                  fSmoothenCounts;        // fill smoothened spectrum with counts
         Bool_t                  fTestMode;              // unfold with unity response for testing
-        Bool_t                  fNoDphi;                // no dphi dependence in unfolding
         Bool_t                  fRawInputProvided;      // input histograms provided, not read from file
         Double_t                fEventPlaneRes;         // event plane resolution for current centrality
         Bool_t                  fUseDetectorResponse;   // add detector response to unfolding
         Bool_t                  fUseDptResponse;        // add dpt response to unfolding
         Bool_t                  fTrainPower;            // don't clear the params of fPower for call to Make
-                                                        // might give more stable results, but also introduces
+                                                        // might give more stable results, but possibly introduces
                                                         // a bias / dependency on previous iterations
-        Bool_t                  fInOutUnfolding;        // do the unfolding in in and out of plane orientation
+        Bool_t                  fDphiUnfolding;         // do the unfolding in in and out of plane orientation
+        Bool_t                  fDphiDptUnfolding;      // do the unfolding in dphi and dpt bins (to fit v2)
+        Bool_t                  fExLJDpt;               // exclude randon cones with leading jet
         // members, set internally
         TProfile*               fRMSSpectrumIn;         // rms of in plane spectra of converged unfoldings
         TProfile*               fRMSSpectrumOut;        // rms of out of plane spectra of converged unfoldings
index f877cb1..5b7a7f7 100644 (file)
-\documentclass[12pt]{book}
+\documentclass[a5paper]{book}
+\usepackage[nottoc,numbib]{tocbibind}
+\usepackage[english]{babel}
+\usepackage[pdftex]{graphics} %enable pdfgraphics support necessary for pdflatex
 \usepackage{graphicx}
-\usepackage[pdftex,bookmarks]{hyperref}
-%-----------------------------------------------------------------------
-% Extra commands
-\include{commands}
-%-----------------------------------------------------------------------
-%\includeonly{chapter4,biblio}
-%-----------------------------------------------------------------------
+\usepackage{tikz}
+\usepackage[small,bf]{caption}
+\usepackage{footmisc}
+\usepackage{sidecap}
+\usepackage{fancyvrb}
+\usepackage{color}
+\usepackage{rotating}
+\makeatletter
+\setlength{\@fptop}{0pt}
+\setlength{\@fpsep}{5pt}
+\makeatother
+\renewcommand{\dotfill}{\leaders\hbox to 5pt{\hss.\hss}\hfill}
+\usepackage[hmargin=1.2cm,vmargin=1.5cm]{geometry} %set custom margins
+\usepackage{fancyhdr} %enable hyperlinking, highlight page numbers
+\usepackage{amsmath}
+\usepackage{amssymb}
+\numberwithin{equation}{subsection}
+\usepackage{amsfonts}
+\usepackage{mathrsfs}
+\usepackage{tabularx}
+\usepackage{lastpage}
+\usepackage{subfigure}
+\usepackage{multicol}
+\usepackage{lipsum}
+\usepackage{xspace}
+\usepackage{latexsym} 
+\usepackage{verbatim}
+\usepackage{listings}
+\usepackage{color}
+\usepackage[usenames,dvipsnames]{xcolor}
+\definecolor{listinggray}{gray}{0.9}
+\definecolor{lbcolor}{rgb}{0.9,0.9,0.9}
+\usepackage{listings}
+\usepackage{color}
+
+\definecolor{mygreen}{rgb}{0,0.6,0}
+\definecolor{mygray}{rgb}{0.5,0.5,0.5}
+\definecolor{mymauve}{rgb}{0.58,0,0.82}
+
+\lstset{ %
+   backgroundcolor=\color{white},   % choose the background color; you must add \usepackage{color} or \usepackage{xcolor}
+   basicstyle=\footnotesize,        % the size of the fonts that are used for the code
+   breakatwhitespace=false,         % sets if automatic breaks should only happen at whitespace
+   breaklines=true,                 % sets automatic line breaking
+%   captionpos=b,                    % sets the caption-position to bottom
+   commentstyle=\color{mygreen},    % comment style
+   deletekeywords={...},            % if you want to delete keywords from the given language
+   escapeinside={\%*}{*)},          % if you want to add LaTeX within your code
+   extendedchars=true,              % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8
+   frame=single,                    % adds a frame around the code
+   keepspaces=true,                 % keeps spaces in text, useful for keeping indentation of code (possibly needs columns=flexible)
+   keywordstyle=\color{blue},       % keyword style
+   language=Octave,                 % the language of the code
+   morekeywords={*,...},            % if you want to add more keywords to the set
+   numbers=left,                    % where to put the line-numbers; possible values are (none, left, right)
+   numbersep=5pt,                   % how far the line-numbers are from the code
+   numberstyle=\tiny\color{mygray}, % the style that is used for the line-numbers
+   rulecolor=\color{black},         % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here))
+   showspaces=false,                % show spaces everywhere adding particular underscores; it overrides 'showstringspaces'
+   showstringspaces=false,          % underline spaces within strings only
+   showtabs=false,                  % show tabs within strings adding particular underscores
+   stepnumber=1,                    % the step between two line-numbers. If it's 1, each line will be numbered
+   stringstyle=\color{mymauve},     % string literal style
+   tabsize=2,                       % sets default tabsize to 2 spaces
+%   title=\lstname                   % show the filename of files included with \lstinputlisting; also try caption instead of title
+}
+% tweak the the style a little bit
+\lstdefinestyle{customc}{
+%   belowcaptionskip=1\baselineskip,
+   breaklines=true,
+   frame=single,
+   xleftmargin=\parindent,
+   language=C++,
+%   showstringspaces=false,
+   basicstyle=\footnotesize\ttfamily,
+   keywordstyle=\bfseries\color{green!40!black},
+   commentstyle=\itshape\color{mygray},
+   identifierstyle=\color{blue},
+   stringstyle=\color{orange},
+   framexleftmargin=5mm,
+}
+
+\lstset{escapechar=@,style=customc}
+\usepackage[sort&compress,numbers]{natbib}
+
+\pagestyle{fancy} %enable custom headeres
+\fancyhead{}
+\fancyfoot{} %clear predefined headers, then create your own:
+\lhead{\footnotesize AliROOT Flow Package manual and documentation}
+\rhead{\footnotesize The FLOW team}
+\rfoot{\footnotesize Page \thepage\ of \pageref{LastPage}}
+\lfoot{\footnotesize \rightmark}
+\usepackage[bookmarks=true,
+            bookmarksopen=true, 
+            bookmarksnumbered=true, pdftex,
+            hyperindex=true,linktocpage]{hyperref}
+\hypersetup{
+    bookmarks=true,         % show bookmarks bar?
+    unicode=false,          % non-Latin characters in Acrobat’s bookmarks
+    pdftoolbar=true,        % show Acrobat’s toolbar?
+    pdfmenubar=true,        % show Acrobat’s menu?
+    pdffitwindow=false,     % window fit to page when opened
+    pdfstartview={FitH},    % fits the width of the page to the window
+    pdftitle={AliROOT Flow Package manual and documentation},    % title
+    pdfauthor={Redmer Alexander Bertens},     % author
+    pdfsubject={AliROOT Flow Package manual and documentation},   % subject of the document
+    pdfcreator={Redmer Alexander Bertens},   % creator of the document
+    pdfproducer={Redmer Alexander Bertens}, % producer of the document
+    pdfkeywords={ALICE} {AliROOT} {ROOT} {manual} {flow} {package} {software}, % list of keywords
+    pdfnewwindow=true,      % links in new window
+    colorlinks=true,       % false: boxed links; true: colored links
+    linkcolor=blue,          % color of internal links
+    citecolor=green,        % color of links to bibliography
+    filecolor=blue,      % color of file links
+    urlcolor=blue           % color of external links
+}
+\usepackage{lineno}
+\usepackage{makeidx}
+\makeindex
+
+\sloppy
+\fancypagestyle{plain}
+ {
+ \fancyhead{}
+ \fancyfoot{}
+ } % clear header and footer of plain page because of ToC
+
+ % tikxstyle for integration of tikz graphics
+\tikzstyle{every node}=[font=\tiny]
+\tikzstyle{decision} = [rectangle, draw, fill=yellow!20, 
+    rounded corners, minimum height=1.5em, minimum width=2.5em, inner sep=3pt]
+\tikzstyle{decision2} = [rectangle, draw, fill=red!20, 
+    rounded corners, minimum height=1.5em, minimum width=2.5em, inner sep=3pt]
+\tikzstyle{decision3} = [rectangle, draw, fill=green!20, 
+    rounded corners, minimum height=1.5em, minimum width=2.5em, inner sep=3pt]
+\tikzstyle{block} = [rectangle, draw, fill=blue!20, 
+    rounded corners, minimum height=1.5em, minimum width=2.5em, inner sep=3pt]
+\tikzstyle{line} = [draw, -latex']
+\tikzstyle{cloud} = [draw, ellipse,fill=red!20, 
+    minimum height=2em]
+\tikzstyle{sdecision} = [rectangle, draw, fill=yellow!20, 
+    rounded corners, minimum height=1.5em, minimum width=2.5em, inner sep=3pt]
+\tikzstyle{sblock} = [rectangle, draw, fill=blue!20, 
+    rounded corners, minimum height=1.5em, minimum width=2.5em, inner sep=3pt]
+\tikzstyle{scloud} = [draw, ellipse,fill=red!20, 
+    minimum height=2em]
+\tikzstyle{ycloud} = [draw, ellipse,fill=yellow!20, 
+    minimum height=2em]
+
+\renewcommand{\thefootnote}{\fnsymbol{footnote}}
+\frontmatter
+
+\linenumbers
 \begin{document}
 
-%\layout
-%\pagenumbering{roman}
-\include{chapters/FrontMatter}
-%-----------------------------------------------------------------------
+\noindent
+\begin{center}
+       \vspace*{1.5cm}
+       {\LARGE \bf The FLOW Analysis Package}\\
+               
+       \vspace{1.5cm}
+       \begin{figure}[hbt]
+               \includegraphics[width=1.\textwidth]{figs/daVinci.png}
+       \end{figure}
+               
+       \vspace{1.5cm}
+       \noindent
+       {\large \bf a short writeup}\\
+       \today\\
+\end{center}
+\vfill
+\noindent
+Various authors, edited by Redmer Alexander Bertens (\texttt{rbertens@cern.ch})
+
+\clearpage
+\thispagestyle{empty}
+\pagenumbering{roman}
 \tableofcontents
+\renewcommand{\thefootnote}{\alph{footnote}}
+\mainmatter
+\chapter{Introduction}
+The intro to everything.
 %-----------------------------------------------------------------------
-%\include{chapters/Introduction}
-%-----------------------------------------------------------------------
-\include{chapters/QuickStart}
-%-----------------------------------------------------------------------
-\include{chapters/FlowEvent}
-%-----------------------------------------------------------------------
-\include{chapters/Program}
-%-----------------------------------------------------------------------
-\include{chapters/Methods}
-%-----------------------------------------------------------------------
-\include{chapters/MonteCarlo}
-%-----------------------------------------------------------------------
-\include{chapters/ScalarProduct}
-%-----------------------------------------------------------------------
-\include{chapters/GeneratingFunctionCumulants}
-%-----------------------------------------------------------------------
-\include{chapters/QVectorCumulants}
-%-----------------------------------------------------------------------
-\include{chapters/LeeYangZeros}
-%-----------------------------------------------------------------------
-\include{chapters/LeeYangZerosEP}
+\chapter{A Quick Start}
+\section{The flow package}
+\label{quickstart}
+The \texttt{ALICE flow package}\index{flow package}\index{ALICE flow package|see{flow package}}\footnote{The ALICE flow package is part of AliROOT\index{AliROOT}, the ALICE extension of the ROOT framework, which can be obtained from \href{http://git.cern.ch/pub/AliRoot}{http://git.cern.ch/pub/AliRoot}. The flow package itself is located in the folder \texttt{\$ALICE\_ROOT/PWG/FLOW/}, where \texttt{\$ALICE\_ROOT} refers to the source directory of AliROOT.} 
+contains most known flow analysis methods.  In this chapter we give a few examples how to setup an
+analysis for the most common cases. The chapters that follow provide more detailed information on the structure of the code 
+and settings of the various flow methods. 
+This write-up is however not a complete listing of the methods, for this the reader is referred to the header files.
+\section{On the fly}
+The macro \texttt{Documentation/examples/runFlowOnTheFlyExample.C}\index{On the fly}\index{runFlowOnTheFlyExample.C}\footnote{In aliroot, this macro can be found at \\ \texttt{\$ALICE\_ROOT/PWGCF/FLOW/Documentation/examples/runFlowOnTheFlyExample}} is a basic example of how the flow package works. 
+In this section we explain the main pieces of that macro.
+\begin{enumerate}
+       \item To use the flow code the flow library needs to be loaded. In\index{libraries, AliROOT} \texttt{AliROOT}:
+       \begin{lstlisting}
+gSystem->Load("libPWGflowBase");\end{lstlisting}
+       In \texttt{root} additional libraries need to be loaded\index{libraries, ROOT}: 
+       \begin{lstlisting}
+gSystem->Load("libGeom");
+gSystem->Load("libVMC");
+gSystem->Load("libXMLIO");
+gSystem->Load("libPhysics");
+gSystem->Load("libPWGflowBase");\end{lstlisting}
+       \item We need to instantiate the flow analysis methods which we want to use. In this example we will
+       instantiate two methods: the first which calculates the flow versus the reaction plane of the Monte Carlo, which is our reference value\index{reference value}\index{Monte Carlo Event Plane} (see section \ref{MC}), and second the so called Q-cumulant method (see section \ref{qvc}).
+       \begin{lstlisting}
+AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
+AliFlowAnalysisWithQCumulants *qc = new AliFlowAnalysisWithQCumulants();\end{lstlisting}
+       \item Each of the methods needs to initialize\index{initialize methods} (e.g. to define the histograms): 
+       \begin{lstlisting}
+mcep->Init(); 
+qc->Init();\end{lstlisting}
+       \item To define the particles we are going to use as Reference Particles\index{reference particles} \index{RP|see{reference particles}} (RP's, particles 
+       used for the {\bf Q} vector) and the Particles Of Interest\index{particles of interest} \index{POI|see{particles of interest}} (POI's, the particles of which 
+       we calculate the differential flow) we have to define two track cut objects\index{track cut object, simple}:
+       \begin{lstlisting}
+AliFlowTrackSimpleCuts *cutsRP = new AliFlowTrackSimpleCuts();
+AliFlowTrackSimpleCuts *cutsPOI = new AliFlowTrackSimpleCuts();
+cutsPOI->SetPtMin(0.2);
+cutsPOI->SetPtMax(2.0);        \end{lstlisting}
+       \item Now we are ready to start the analysis. For a quick start we make an event on the fly, tag the reference particles and particles of interest  and pass it to the two flow methods. \\
+       \begin{lstlisting}
+for(Int_t i=0; i<nEvts; i++) { 
+      // make an event with mult particles 
+      AliFlowEventSimple* event = new AliFlowEventSimple(mult,AliFlowEventSimple::kGenerate);
+      // modify the tracks adding the flow value v2
+      event->AddV2(v2);
+      // select the particles for the reference flow
+      event->TagRP(cutsRP);
+      // select the particles for differential flow
+      event->TagPOI(cutsPOI);
+      // do flow analysis with various methods:
+      mcep->Make(event);
+      qc->Make(event);
+}\end{lstlisting}
+       \item To fill the histograms which contain the final results we have to call Finish for each method:
+       \begin{lstlisting}
+mcep->Finish(); 
+qc->Finish(); \end{lstlisting}
+       \item This concludes the analysis and now we can write the results into a file:
+       \begin{lstlisting}
+TFile *outputFile = new TFile("AnalysisResults.root","RECREATE");
+mcep->WriteHistograms();
+qc->WriteHistograms();\end{lstlisting}
+\end{enumerate}
+
+\section{What is in the output file?}
+\index{output file}Now we have written the results into a file, but what is in there?
+
+\section{Reading events from file}
+The macro \texttt{Documentation/examples/runFlowReaderExample.C} \index{runFlowReaderExample.C}is an example how to setup a flow analysis if the events are already generated and
+for example are stored in ntuples.
+\section{A simple flow analysis in ALICE using Tasks}
+The macro \texttt{Documentation/examples/runFlowTaskExample.C} \index{runFlowTaskExample.C}is an example how to setup a flow analysis using the full ALICE Analysis Framework.
+
 %-----------------------------------------------------------------------
-\include{chapters/FitQDistribution}
+
+\chapter{The Flow Event}
+\label{flowevent}
+Here we describe the flowevent, flowtracks, general cuts and cuts for RPs POIs.
+OntheFly, AfterBurner. Filling with ESD, AOD, Ntuples, etc.  
 %-----------------------------------------------------------------------
-\include{chapters/Summary}
+
+\chapter{The Program}
+\label{The Program}
+\begin{center}
+       \includegraphics[width=1.\textwidth]{figs/flowChart.png}
+\end{center}
+Here we describe the program.
 %-----------------------------------------------------------------------
-\include{chapters/bibliography}
+
+\chapter{Methods}
+\label{Methods}
+The flow package aims at providing the user with most of the known flow analysis methods. Detailed theoretical overview of the methods can be found in the following papers, which are included in the folder \texttt{\$ALICE\_ROOT/PWGCF/FLOW/Documentation/otherdocs/}
+\begin{itemize}
+       \item The Monte-Carlo Truth
+       \item Scalar Product Method \\ \hspace*{1cm} \texttt{EventPlaneMethod/FlowMethodsPV.pdf}
+       \item Generating Function Cumulants \\ \hspace*{1cm} \texttt{GFCumulants/Borghini\_GFCumulants\_PracticalGuide.pdf}
+       \item Q-vector Cumulant method \\ \hspace*{1cm} \texttt{QCumulants/QCpaperdraft.pdf} 
+        \item Lee-Yang Zero Method \\ \hspace*{1cm} \texttt{LeeYangZeroes/Borghini\_LYZ\_PracticalGuide.pdf}
+       \item Lee-Yang Zero Method \\ \hspace*{1cm} \texttt{LeeYangZeroesEP/LYZ\_RP.pdf}
+\end{itemize}
+The structure of this  chapter is as follows: of each of the available methods a short description is given in the \texttt{theory} subsection (for more detailed information, see the papers listed above) followed by details which are specific to the implementation in the subsection \texttt{implementation}. Caveats, possible issues, etc, are listed in the \texttt{caveats} subsections. 
 %-----------------------------------------------------------------------
-\include{chapters/appendix1}
+
+\section{The Monte-Carlo Truth}
+\label{MC}
+Here we describe the implementation of the monte-carlo truth.
 %-----------------------------------------------------------------------
-\end{document}
 
+\section{Scalar Product Method}
+\label{SP}
+\subsection{Theory}
+\subsubsection{The scalar product method}
+The scalar product method - as well as the Q-cumulant method which will be described later - does not depend on explicit construction of an (sub)event plane, but estimates $v_n$ directly from multi-particle correlations. 
+
+To do so, firstly all particles in an event are labeled either as \emph{reference particles} (RP's) or \emph{particles of interest} (POI's). The RP and POI selections are in turn divided into sub-events, which are again taken from different $\eta$ ranges, in analogy to the approach taken for the event plane method. Each POI is correlated with a sub-event Q-vector from the RP selection, which allows for the calculation of $v_n$ without any explicit reconstruction of an event plane\cite{sp}. 
 
+The reason for the division into RP's and POI's is the fact that the two particle correlator of POI's,
+\begin{equation}\label{unstable}
+       v_n^{POI} = \sqrt{ \left< e^{i n (\phi^{POI}_i - \phi^{POI}_j)} \right> }
+\end{equation}
+is generally not stable statistically. Introducing reference flow, \ref{unstable} can be rewritten as
+\begin{equation}\label{stable}
+       v_n^{POI} = \frac{ \left< e^{i n (\phi^{POI}_i - \phi^{RP}_j)} \right>}{\left< e^{i n (\phi^{RP}_i - \phi^{RP}_j)} \right>} = \frac{v_n^{POI} v_n^{RP} }{\sqrt{v_n^{RP} v_n^{RP}}}.
+\end{equation}
+By taking an abundant particle source as RP's - in the case of this study the RP selection comprises all charged particles - both correlators in \ref{stable} are statistically stable. 
 
+\paragraph{The scalar product method}
+In the scalar product method, POI's $u_k$,
+\begin{equation}\label{spderiv}
+       u_k = e^{i n \phi_k},
+\end{equation}
+are correlated with $Q^*_a$, the complex-conjugate Q-vector built from RP's in a given sub-event $a$. First, the scalar product of $u_k$ and $Q^*_a$ is taken,
+\begin{equation}
+       u_k \cdotp \sum_{\substack{ j=1,\\j \neq k}}^{M_{RP, a}} u^*_j
+\end{equation}
+where $M_{RP, a}$ denotes RP multiplicity for a given sub-event $a$ and the inequality $j \neq k$ removes auto-correlations. From this, differential $v_n$ of POI's ($v_n^{\prime}$) and $v_n$ of RP's ($v_n^a$) in sub-event $a$ can be obtained in a straightforward way from the correlation of POI's and RP's:
+\begin{equation}\label{rppoisp}
+       \left< u \cdotp Q^*_a \right> = \frac{1}{M_{RP, a}-k}\sum_{i=k}^{M_{RP, a}} \left( u_k \sum_{\substack{ j=1,\\j \neq k}}^{M_{RP, a}} u^*_j \right)
+\end{equation}
+where POI multiplicity is expressed in terms of $M_{RP, a}$; $M_{POI} = M_{RP, a} - k$. Since for any function $f(x)$ and constant $a$
+\begin{equation}
+       \sum a f(x) = a \sum f(x) 
+\end{equation}
+\ref{rppoisp} can be rewritten as
+\begin{align}\label{sp_part2}
+       \left< u \cdotp Q^*_a \right> & = \frac{1}{M_{RP, a}-k}\sum_{i=k}^{M_{RP, a}} e^{i n [\phi_k - \Psi_n]} \sum_{j=1}^{M_{RP, a}} e^{- i n [\phi_j - \Psi_n]} \\
+                                     & = M_{RP, a} v^{\prime}_n v_n^a \nonumber                                                                                   
+\end{align}
+where in the last step of \ref{sp_part2} it has been used that
+\begin{equation}
+       v_n = \frac{\ds \sum_i^M e^{in[\phi_i - \Psi_n]} }{M}.
+\end{equation}
 
+To obtain the estimate of $v_n$, one must still disentangle the reference flow contribution from the event averaged correlation given in \ref{rppoisp}. 
+Proceeding in a fashion similar to that presented in equation \ref{rppoisp}, it can be shown that
+\begin{equation}
+       \left< \frac{Q_a}{M_a}\cdotp\frac{Q^*_b}{M_b} \right> = \left< v_n^{a}v_n^{b} \right>
+\end{equation}
+where $Q_a, Q_b$ are the Q-vectors of RP's in sub-event $a, b$. Under the assumption that
+\begin{equation}\label{zerononflow}
+       \left< v_n^2 \right> = \left< v_n \right>^2,
+\end{equation}
+- an assumption which will be spoiled in the case of flow fluctuations - and requiring that the $v_n$ estimates in both sub-events are equal, one simply evaluates
+\begin{equation}\label{sp}
+       v_n^{\prime} = \frac{\left< \left< u \cdotp \frac{Q^*_a}{M_a} \right>\right>}{\sqrt{ \left< \frac{Q_a}{M_a}\cdotp\frac{Q^*_b}{M_b} \right>}}
+\end{equation}
+to obtain $v_n^{a}$. For equal multiplicity sub-events $M_a = M_b$, \ref{sp} is simplified to
+\begin{equation}\label{spa}
+       v_n^{\prime} = \frac{\left< \left< u \cdotp Q^*_a \right>_t \right>}{ \sqrt{ \left< Q_a\cdotp Q^*_b \right>}}.
+\end{equation}
+$v_n^b$ can be obtained by switching indices $a$ and $b$ in expressions \ref{sp} and \ref{spa}, and should equal $v_n^a$.
+This principle can be generalized straightforwardly to allow for a selection of RP's which has been divided into three subevents. 
+\begin{align}
+       v_n^a & = \frac{ \left< \left< u \cdotp \frac{Q^*_a}{M_a} \right> \right>}{\sqrt{\left< v_n^{\prime a} v_n^{\prime b} \right> \left< v_n^{\prime a}v_n^{\prime c} \right> / \left< v_n^{\prime b} v_n^{\prime c} \right>}}                                                    \\
+             & = \frac{\left< \left< u \cdotp \frac{Q^*_a}{M_a} \right> \right>}{\sqrt{\left< \frac{Q_a}{M_a} \cdotp \frac{Q_{b}^{*}}{M_b} \right> \left< \frac{Q_a}{M_a} \cdotp \frac{Q_{c}^{*}}{M_c} \right> / \left< \frac{Q_b}{M_b} \cdotp \frac{Q_c^*}{M_c}} \right>} \nonumber 
+\end{align}
+where cyclic permutation of $a, b, c$ (in analogy to the switching of indices in \ref{3evsp} gives the estimates of $v_n^b$ and $v_n^c$.\textcolor{red}{[insert some discussion here: is this result actually true, and some light on va, vb, (vc)]}
+
+\subsection{Implementation}
+
+\subsubsection{Extension to Event Plane method}
+As explained earlier, the event plane analysis results in this study are actually obtained by normalizing the Q-vectors in the scalar product by their length $\vert Q_n \vert$. Consider the following:
+\begin{equation}\label{phase}
+       \frac{Q_n^*}{\vert Q_n^* \vert} = \frac{\vert Q_n^* \vert e^{- i n \Psi_{Q_n}}}{\vert Q_n^* \vert} = e^{- i n \Psi_{Q_n}}.
+\end{equation}
+For a full event, the enumerator of \ref{sp} can be expressed as
+\begin{equation}\label{epsp1}
+       \left< \left< u \cdotp e^{- i n \Psi_{Q_n}} \right> \right> = \left< \left< e^{i n \phi_i} \cdotp e^{- i n \Psi_{Q_n}} \right> \right> \nonumber = \left< \left< e^{i n (\phi_i - \Psi_{Q_n})} \right> \right> = \left< \left< \cos(n [\phi_i - \Psi_{Q_n}]) \right> \right>
+\end{equation}
+which corresponds to the all-event average of \ref{ep_first}. As shown in the previous subsection this expression equals $v_n^{obs}$. 
+
+For normalized Q-vectors, the denominator of \ref{sp} reads (using \ref{phase}):
+\begin{equation}\label{epsp2}
+       \sqrt{\left< \frac{Q_a}{\vert Q_a\vert} \cdotp \frac{Q_b^*}{\vert Q_b^* \vert \right>}} = \sqrt{\left< e^{in[\Psi_{Q_{n_a}} - \Psi_{Q_{n_b}}}\right>} = \sqrt{\left< \cos(n[\Psi_{Q_{n_a}} - \Psi_{Q_{n_b}}]\right>)}
+       \end{equation}
+       from which the event plane resolution can be calculated using \ref{emeancos} or \ref{halfeventresolution}.
+               
+       \subsubsection{Caveats}
+               
+       %-----------------------------------------------------------------------
+               
+       \section{Generating Function Cumulant Method}
+       \label{GFC}
+       Here we describe the generating function cumulant method and how it is implemented.
+       %-----------------------------------------------------------------------
+               
+       \section{Q-vector Cumulant Method}
+       \label{qvc}
+       \subsection{Theory}
+       The Q-cumulant (QC) method\footnote{The overview given in this section is inspired by \cite{bilandzic-2011-83}, for further reading the reader is referred there. A full derivation of results that are relevant in this study is given in appendix \ref{qcapp}.} uses multi-particle correlations to estimate $v_n$ estimates for RP's and POI's, but does not limit itself to two-particle correlations. Although higher-order Q-cumulant calculations are available, this section will discuss the method using two- and four-particle correlations. 
+               
+       %\paragraph{Multiparticle correlations and anisotropic flow}
+       Multi-particle correlations in the QC method are expressed in terms of cumulants, which are the the expectation values of correlation terms in joint probability density functions. Consider the following: if two observables $f$ for particles $x_i$ and $x_j$ are correlated, the joint probability $f(x_i, x_j)$ is the sum of the factorization of the constituent probabilities and a covariance term:
+       \begin{equation}
+               f(x_i, x_j) = f(x_i)f(x_j) + f_c (x_i, x_j)
+       \end{equation}
+       When taking as an observable azimuthal dependence,
+       \begin{align}
+               x_i \equiv e^{i n \phi_i}, \hspace{.5in} x_j \equiv e^{i n \phi_j} 
+       \end{align}
+       the two-particle cumulant is expressed as the covariance of the expectation value:
+       \begin{equation}
+               E_C(e^{i n[\phi_i - \phi_j]}) = E(e^{[i n(\phi_i - \phi_j]}) - E(e^{in[\phi_i]})E(e^{in[- \phi_j]}).
+       \end{equation}
+       Symmetry arguments (along the lines of those given in appendix \ref{sec:harmonic_derivation}) dictate that the product of separate expectation values is equals zero, from which a familiar expression for the two-particle correlation is obtained:
+       \begin{align}
+               E_C(e^{i n[\phi_i - \phi_j]}) & = E(e^{in[\phi_i]})E(e^{in[- \phi_j]})                                   \\
+                                             & = \left< e^{in[\phi_i]}\right> \left< e^{in[- \phi_j]} \right> \nonumber \\
+                                             & = \left< e^{in[\phi_i - \phi_j]} \right> \nonumber                       \\
+                                             & = \left< 2 \right>, \nonumber                                            
+       \end{align}
+       the all-event average of which is denoted by
+       \begin{equation}\label{twocul}
+               c_n\{2\} = \left< \left< 2 \right> \right>
+       \end{equation}
+       where $c_n\{2\}$ is called the two-particle cumulant. For the four-particle case, one proceeds likewise:
+       \begin{align}
+               E_c(e^{in[\phi_i + \phi_j - \phi_k -\phi_l]}) & = E(e^{in[\phi_i + \phi_j - \phi_k -\phi_l]})                     \\
+                                                             & - E(e^{in[\phi_i - \phi_k]})E(e^{in[\phi_j - \phi_l]})\nonumber   \\
+                                                             & - E(e^{in[\phi_i - \phi_l]})E(e^{in[\phi_j - \phi_k]}). \nonumber 
+       \end{align}
+       The four-particle cumulant can be expressed in terms of two- and four-particle correlations as well,
+       \begin{equation}\label{fourcul}
+               c_n\{4\} = \left< \left< 4 \right> \right> - 2 \left< \left< 2 \right> \right>^2.
+       \end{equation}
+       From \ref{twocul} and \ref{fourcul} it follows that $v_n$ harmonics are related to cumulants following
+       \begin{align}\label{refFlowFromCumulants2nd}
+               v_n\{2\} & = \sqrt{c_n\{2\}}                \\
+               v_n\{4\} & = \sqrt[4]{-c_n\{4\}} \nonumber. 
+       \end{align}
+       where $v_n\{2\}$, $v_n\{4\}$ denote flow estimates obtained from two- and four-particle correlations.
+               
+       In a fashion similar to that explained in the previous subsection, the Q-cumulant method uses reference flow to obtain a statistically stable estimate of the differential flow of POI's. Differential POI flow, for the two- and four-particle case, can be expressed as
+       \begin{align}
+               d_n\{2\} & = \left< \left< 2^{\prime} \right> \right>                                                                                            \\
+               d_n\{4\} & = \left< \left< 4^{\prime} \right> \right> - 2\cdotp \left< \left< 2^{\prime} \right> \right>\left< \left< 2 \right> \right>\nonumber 
+       \end{align}
+       where $d_n\{2\}, d_n\{4\}$ denotes the two-, four-particle differential flow and the $^{\prime}$ is used as an indicator for differential (p$_t$ dependent) results. Disentangling from this the reference flow contributions, one is left with the final expression for the estimate of differential $v_n$ for POI's:
+       \begin{align}
+               v_n^{\prime}\{2\} & = \frac{d_n\{2\}}{\sqrt{c_n\{2\}}}               \\
+               v_n^{\prime}\{4\} & = - \frac{d_n\{4\}}{(-c_n\{2\})^{3/4}}.\nonumber 
+       \end{align}
+               
+       \subsection{Implementation}
+               
+       Here we describe the Q-vector cumulant method and how it is implemented.
+       %-----------------------------------------------------------------------
+               
+       \section{Lee-Yang Zero Method}
+       \label{LYZ}
+       Here we describe the Lee-Yang Zero method and how it is implemented.
+       %-----------------------------------------------------------------------
+               
+       \section{Lee-Yang Zero Method}
+       \label{LYZ}
+       Here we describe the Lee-Yang Zero method and how it is implemented.
+       %-----------------------------------------------------------------------
+               
+       \section{Fitting the Q-vector Distribution}
+       \label{qfit}
+       Here we describe how the fitting of the Q-vector distribution is implemented.
+       %-----------------------------------------------------------------------
+               
+       \chapter{Summary}
+       \label{Summary}
+       This sums it all up.
+       %-----------------------------------------------------------------------
+               
+       \begin{thebibliography}{99}
+                               
+                                 
+               \bibitem{Ollitrault:1992bk}
+               J.~Y.~Ollitrault,
+               %``Anisotropy as a signature of transverse collective flow,''
+               Phys.\ Rev.\ D {\bf 46} (1992) 229.
+               %%CITATION = PHRVA,D46,229;%%
+                               
+               \bibitem{Danielewicz:1999vh}
+               P.~Danielewicz,
+               %``Flow and equation of state in heavy-ion collisions,''
+               Nucl.\ Phys.\ A {\bf 661} (1999) 82.
+               %  [arXiv:nucl-th/9907098].
+               %%CITATION = NUCL-TH 9907098;%%
+                               
+               \bibitem{Rischke:1996nq}
+               D.~H.~Rischke,
+               %``Hydrodynamics and collective behavior in relativistic nuclear
+               %collisions,''
+               Nucl.\ Phys.\ A {\bf 610} (1996) 88C.
+               %  [arXiv:nucl-th/9608024].
+               %%CITATION = NUCL-TH 9608024;%%
+                               
+               \bibitem{Ollitrault:1997vz}
+               J.~Y.~Ollitrault,
+               %``Flow systematics from SIS to SPS energies,''
+               Nucl.\ Phys.\ A {\bf 638} (1998) 195.
+               %  [arXiv:nucl-ex/9802005].
+               %%CITATION = NUCL-EX 9802005;%%
+                               
+               \bibitem{Voloshin:1994mz}
+               S.~Voloshin and Y.~Zhang,
+               %``Flow study in relativistic nuclear collisions by Fourier expansion of
+               %Azimuthal particle distributions,''
+               Z.\ Phys.\ C {\bf 70} (1996) 665.
+               %  [arXiv:hep-ph/9407282].
+               %%CITATION = HEP-PH 9407282;%%
+                               
+               %\cite{Ackermann:2000tr}
+               \bibitem{Ackermann:2000tr}  
+               K.~H.~Ackermann {\it et al.}  [STAR Collaboration],
+               %``Elliptic flow in Au + Au collisions at s(N N)**(1/2) = 130-GeV,''
+               Phys.\ Rev.\ Lett.\  {\bf 86} (2001) 402
+               %[arXiv:nucl-ex/0009011].  
+               %%CITATION = NUCL-EX 0009011;%%
+                                 
+               %\cite{Adler:2001nb}
+               \bibitem{Adler:2001nb}  
+               C.~Adler {\it et al.}  [STAR Collaboration],
+               % ``Identified particle elliptic flow in Au + Au collisions at  
+               %s(NN)**(1/2) =  130-GeV,''
+               Phys.\ Rev.\ Lett.\  {\bf 87} (2001) 182301  
+               %[arXiv:nucl-ex/0107003].  
+               %%CITATION = NUCL-EX 0107003;%%
+                               
+               %%% theory RHIC discoveries %%%
+               %\cite{Gyulassy:2004zy}
+               \bibitem{RHIC_Discoveries}
+               T.D.~Lee {\it et al.}, 
+               New Discoveries at RHIC: Case for the Strongly Interacting 
+               Quark-Gluon Plasma. 
+               Contributions from the RBRC Workshop held May 14-15, 2004.
+               Nucl.\ Phys.\ A {\bf 750} (2005) 1-171
+               %%% end theory white-papers
+                               
+                                   
+       \end{thebibliography}
+               
+       %-----------------------------------------------------------------------
+       \chapter*{Appendix I}
+       \label{appendix1}
+       Here we put short pieces of code.
+       %-----------------------------------------------------------------------
+       \printindex
+\end{document}
diff --git a/PWGCF/FLOW/Documentation/chapters/FitQDistribution.tex b/PWGCF/FLOW/Documentation/chapters/FitQDistribution.tex
deleted file mode 100644 (file)
index 8d03aa4..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-\section{Fitting the Q-vector Distribution}
-\label{qfit}
-Here we describe how the fitting of the Q-vector distribution is implemented.
diff --git a/PWGCF/FLOW/Documentation/chapters/FlowEvent.tex b/PWGCF/FLOW/Documentation/chapters/FlowEvent.tex
deleted file mode 100644 (file)
index d05ef69..0000000
+++ /dev/null
@@ -1,4 +0,0 @@
-\chapter{The Flow Event}
-\label{flowevent}
-Here we describe the flowevent, flowtracks, general cuts and cuts for RPs POIs.
-OntheFly, AfterBurner. Filling with ESD, AOD, Ntuples, etc.  
diff --git a/PWGCF/FLOW/Documentation/chapters/FrontMatter.tex b/PWGCF/FLOW/Documentation/chapters/FrontMatter.tex
deleted file mode 100644 (file)
index 58b8992..0000000
+++ /dev/null
@@ -1,20 +0,0 @@
-\thispagestyle{empty}
-\noindent
-\begin{center}
-%\vspace{1.5cm}
-{\LARGE \bf The FLOW Analysis Package}\\
-
-\vspace{1.5cm}
-
-\begin{figure}[hbt]
-    \includegraphics[width=1.\textwidth]{figs/daVinci.png}
-\end{figure}
-
-\vspace{1.5cm}
-\noindent
-{\large \bf a short writeup}\\
-25-10-2010\\
-\end{center}
-
-\clearpage
-\thispagestyle{empty}
diff --git a/PWGCF/FLOW/Documentation/chapters/GeneratingFunctionCumulants.tex b/PWGCF/FLOW/Documentation/chapters/GeneratingFunctionCumulants.tex
deleted file mode 100644 (file)
index edc1982..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-\section{Generating Function Cumulant Method}
-\label{GFC}
-Here we describe the generating function cumulant method and how it is implemented.
diff --git a/PWGCF/FLOW/Documentation/chapters/Introduction.tex b/PWGCF/FLOW/Documentation/chapters/Introduction.tex
deleted file mode 100644 (file)
index 70d0cec..0000000
+++ /dev/null
@@ -1,2 +0,0 @@
-\chapter{Introduction}
-The intro to everything.
\ No newline at end of file
diff --git a/PWGCF/FLOW/Documentation/chapters/LeeYangZeros.tex b/PWGCF/FLOW/Documentation/chapters/LeeYangZeros.tex
deleted file mode 100644 (file)
index 6cc601f..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-\section{Lee-Yang Zero Method}
-\label{LYZ}
-Here we describe the Lee-Yang Zero method and how it is implemented.
diff --git a/PWGCF/FLOW/Documentation/chapters/LeeYangZerosEP.tex b/PWGCF/FLOW/Documentation/chapters/LeeYangZerosEP.tex
deleted file mode 100644 (file)
index 0147b63..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-\section{Lee-Yang Zero Event Plane Method}
-\label{LYZEP}
-Here we describe the Lee-Yang Zero Event Plane method and how it is implemented.
diff --git a/PWGCF/FLOW/Documentation/chapters/Methods.tex b/PWGCF/FLOW/Documentation/chapters/Methods.tex
deleted file mode 100644 (file)
index bc6efcc..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-\chapter{Methods}
-\label{Methods}
-Here we put an intro for  the various methods.
diff --git a/PWGCF/FLOW/Documentation/chapters/MonteCarlo.tex b/PWGCF/FLOW/Documentation/chapters/MonteCarlo.tex
deleted file mode 100644 (file)
index d09f5b3..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-\section{The Monte-Carlo Truth}
-\label{MC}
-Here we describe the implementation of the monte-carlo truth.
\ No newline at end of file
diff --git a/PWGCF/FLOW/Documentation/chapters/Program.tex b/PWGCF/FLOW/Documentation/chapters/Program.tex
deleted file mode 100644 (file)
index 6f1c1f5..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-\chapter{The Program}
-\label{The Program}
-Here we describe the program.
diff --git a/PWGCF/FLOW/Documentation/chapters/QVectorCumulants.tex b/PWGCF/FLOW/Documentation/chapters/QVectorCumulants.tex
deleted file mode 100644 (file)
index 152259e..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-\section{Q-vector Cumulant Method}
-\label{qvc}
-Here we describe the Q-vector cumulant method and how it is implemented.
diff --git a/PWGCF/FLOW/Documentation/chapters/QuickStart.tex b/PWGCF/FLOW/Documentation/chapters/QuickStart.tex
deleted file mode 100644 (file)
index 725ed54..0000000
+++ /dev/null
@@ -1,78 +0,0 @@
-\chapter{A Quick Start}
-\label{quickstart}
-The \textit{ALICE flow package}\footnote{\texttt{http://alisoft.cern.ch/viewvc/trunk/PWG2/FLOW/?root=AliRoot .}} 
-contains most known flow analysis methods.  In this chapter we give a few examples how to setup an
-analysis for the most common cases. The chapters that follow provide more detailed information on the structure of the code 
-and settings of the various flow methods. 
-This write-up is however not a complete listing of the methods, for this the reader is referred to the header files.
-\section{On the fly}
-The macro \texttt{Documentation/examples/runFlowOnTheFlyExample.C} 
- is a basic example of how the flow package works. 
-In this section we explain the main pieces of that macro.
-
-\begin{enumerate}
-\item
-To use the flow code the flow library needs to be loaded. In  \textit{AliRoot}:\\
-\texttt{gSystem->Load("libPWGflowBase");}\\
-In  \textit{root} additional libraries need to be loaded: \\
-\texttt{gSystem->Load("libGeom");}\\
-\texttt{gSystem->Load("libVMC");}\\
-\texttt{gSystem->Load("libXMLIO");}\\
-\texttt{gSystem->Load("libPhysics");}\\
-\texttt{gSystem->Load("libPWGflowBase");}\\
-\item
-We need to instantiate the flow analysis methods which we want to use. In this example we will
-instantiate two methods: the first  which calculates the flow versus the reaction plane of the Monte Carlo, which is our reference value (see section \ref{MC}), 
-and second the so called Q-cumulant method (see section \ref{qvc}).
-\texttt{AliFlowAnalysisWithMCEventPlane *mcep} \\
-\texttt{= new AliFlowAnalysisWithMCEventPlane();}\\
-\texttt{AliFlowAnalysisWithQCumulants *qc}\\
- \texttt{ = new AliFlowAnalysisWithQCumulants();}\\
- \item
- Each of the methods needs to initialize (e.g. to define the histograms): \\
- \texttt{mcep->Init(); }
-\texttt{qc->Init();}\\
-\item
-To define the particles we are going to use as Reference Particles (RP's, particles 
-used for the {\bf Q} vector) and the Particles Of Interest (POI's, the particles of which 
-we calculate the differential flow) we have to define two trackcut objects:\\
-\texttt{AliFlowTrackSimpleCuts *cutsRP = new AliFlowTrackSimpleCuts();}\\
-\texttt{AliFlowTrackSimpleCuts *cutsPOI = new AliFlowTrackSimpleCuts();}\\
-\texttt{cutsPOI->SetPtMin(0.2);}\\
-\texttt{cutsPOI->SetPtMax(2.0);}\\
-\item
-Now we are ready to start the analysis.  
-For a quick start we make an event on the fly, tag the reference particles and particles of interest  and pass it to the two flow methods. \\
-\texttt{for(Int\textunderscore t i=0; i<nEvts; i++) \{}\\
-\texttt{      // make an event with mult particles }\\
-\texttt{      AliFlowEventSimple* event = new AliFlowEventSimple(mult,AliFlowEventSimple::kGenerate);}\\
-\texttt{      // modify the tracks adding the flow value v2}\\
-\texttt{       event->AddV2(v2);}\\
-\texttt{      // select the particles for the reference flow}\\
-\texttt{      event->TagRP(cutsRP);}\\
-\texttt{      // select the particles for differential flow}\\
-\texttt{      event->TagPOI(cutsPOI);}\\
-\texttt{      // do flow analysis with various methods:}\\
-\texttt{      mcep->Make(event);}\\
-\texttt{      qc->Make(event);}\\
-\texttt{    \} // end of for(Int\textunderscore t i=0;i<nEvts;i++)}\\
-\item
-To fill the histograms which contain the final results we have to call Finish for each method:\\
-\texttt{ mcep->Finish(); }  \texttt{ qc->Finish(); }\\
-\item
-This concludes the analysis and now we can write the results into a file:\\
-\texttt{ TFile *outputFile = new TFile("AnalysisResults.root","RECREATE");}\\
-\texttt{ mcep->WriteHistograms();}\\
-\texttt{ qc->WriteHistograms();}\\
-
-\section{What is in the output file?}
-Now we have written the results into a file, but what is in there?
-
-\section{Reading events from file}
-The macro \texttt{Documentation/examples/runFlowReaderExample.C} is an example how to setup a flow analysis if the events are already generated and
-for example are stored in ntuples.
-\section{A simple flow analysis in ALICE using Tasks}
-The macro \texttt{Documentation/examples/runFlowTaskExample.C} is an example how to setup a flow analysis using the full ALICE Analysis Framework.
-\end{enumerate}
diff --git a/PWGCF/FLOW/Documentation/chapters/ScalarProduct.tex b/PWGCF/FLOW/Documentation/chapters/ScalarProduct.tex
deleted file mode 100644 (file)
index ac61108..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-\section{Scalar Product Method}
-\label{SP}
-Here we describe the scalar product method and how it is implemented.
diff --git a/PWGCF/FLOW/Documentation/chapters/Summary.tex b/PWGCF/FLOW/Documentation/chapters/Summary.tex
deleted file mode 100644 (file)
index 6b165c3..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-\chapter{Summary}
-\label{Summary}
-This sums it all up.
\ No newline at end of file
diff --git a/PWGCF/FLOW/Documentation/chapters/appendix1.tex b/PWGCF/FLOW/Documentation/chapters/appendix1.tex
deleted file mode 100644 (file)
index 94c4d3a..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-\chapter*{Appendix I}
-\label{appendix1}
-Here we put short pieces of code.
diff --git a/PWGCF/FLOW/Documentation/chapters/bibliography.tex b/PWGCF/FLOW/Documentation/chapters/bibliography.tex
deleted file mode 100644 (file)
index 3fbad55..0000000
+++ /dev/null
@@ -1,69 +0,0 @@
-\begin{thebibliography}{99}
-
-  
-\bibitem{Ollitrault:1992bk}
-  J.~Y.~Ollitrault,
-  %``Anisotropy as a signature of transverse collective flow,''
-  Phys.\ Rev.\ D {\bf 46} (1992) 229.
-  %%CITATION = PHRVA,D46,229;%%
-
-\bibitem{Danielewicz:1999vh}
-  P.~Danielewicz,
-  %``Flow and equation of state in heavy-ion collisions,''
-  Nucl.\ Phys.\ A {\bf 661} (1999) 82.
-  %  [arXiv:nucl-th/9907098].
-  %%CITATION = NUCL-TH 9907098;%%
-
-\bibitem{Rischke:1996nq}
-  D.~H.~Rischke,
-  %``Hydrodynamics and collective behavior in relativistic nuclear
-  %collisions,''
-  Nucl.\ Phys.\ A {\bf 610} (1996) 88C.
-  %  [arXiv:nucl-th/9608024].
-  %%CITATION = NUCL-TH 9608024;%%
-
-\bibitem{Ollitrault:1997vz}
-  J.~Y.~Ollitrault,
-  %``Flow systematics from SIS to SPS energies,''
-  Nucl.\ Phys.\ A {\bf 638} (1998) 195.
-%  [arXiv:nucl-ex/9802005].
-  %%CITATION = NUCL-EX 9802005;%%
-
-\bibitem{Voloshin:1994mz}
-  S.~Voloshin and Y.~Zhang,
-  %``Flow study in relativistic nuclear collisions by Fourier expansion of
-  %Azimuthal particle distributions,''
-  Z.\ Phys.\ C {\bf 70} (1996) 665.
-  %  [arXiv:hep-ph/9407282].
-  %%CITATION = HEP-PH 9407282;%%
-
-  %\cite{Ackermann:2000tr}
-\bibitem{Ackermann:2000tr}  
-  K.~H.~Ackermann {\it et al.}  [STAR Collaboration],
-  %``Elliptic flow in Au + Au collisions at s(N N)**(1/2) = 130-GeV,''
-  Phys.\ Rev.\ Lett.\  {\bf 86} (2001) 402
-  %[arXiv:nucl-ex/0009011].  
-  %%CITATION = NUCL-EX 0009011;%%
-  
-    %\cite{Adler:2001nb}
-\bibitem{Adler:2001nb}  
-  C.~Adler {\it et al.}  [STAR Collaboration],
-  % ``Identified particle elliptic flow in Au + Au collisions at  
-  %s(NN)**(1/2) =  130-GeV,''
-  Phys.\ Rev.\ Lett.\  {\bf 87} (2001) 182301  
-  %[arXiv:nucl-ex/0107003].  
-  %%CITATION = NUCL-EX 0107003;%%
-
- %%% theory RHIC discoveries %%%
-  %\cite{Gyulassy:2004zy}
-\bibitem{RHIC_Discoveries}
-  T.D.~Lee {\it et al.}, 
-  New Discoveries at RHIC: Case for the Strongly Interacting 
-  Quark-Gluon Plasma. 
-  Contributions from the RBRC Workshop held May 14-15, 2004.
-  Nucl.\ Phys.\ A {\bf 750} (2005) 1-171
-  %%% end theory white-papers
-
-    
-\end{thebibliography}
-
diff --git a/PWGCF/FLOW/Documentation/figs/flowChart.png b/PWGCF/FLOW/Documentation/figs/flowChart.png
new file mode 100644 (file)
index 0000000..bbc51ba
Binary files /dev/null and b/PWGCF/FLOW/Documentation/figs/flowChart.png differ