PWG/FLOW/Base/AliFlowAnalysisWithMSP.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 #include "AliFlowAnalysisWithMSP.h"
  17
  18 #include "AliFlowVector.h"
  19 #include "AliFlowMSPHistograms.h"
  20 #include "AliFlowEventSimple.h"
  21 #include "AliFlowTrackSimple.h"
  22 #include "AliFlowCommonConstants.h"
  23 #include "AliFlowCommonHist.h"
  24 #include "AliFlowCommonHistResults.h"
  25
  26 #include <TList.h>
  27 #include <TH1D.h>
  28 #include <TProfile.h>
  29 #include <TVector2.h>
  30 #include <TDirectoryFile.h>
  31 #include <TMath.h>
  32
  33 #include <iostream>
  34 #include <iomanip>
  35
  36 AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP()
  37    : TNamed(), fHarmonic(2), fNUA(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0),
  38    fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0),
  39    fPtStatistics(0), fEtaStatistics(0),
  40    fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0)
  41 {
  42    // Default constructor. Intended for root IO purposes only
  43 }
  44
  45 AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(TDirectoryFile *file)
  46    : TNamed(), fHarmonic(2), fNUA(kFALSE), fUseCommonConstants(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0),
  47    fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0),
  48    fPtStatistics(0), fEtaStatistics(0),
  49    fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0)
  50 {
  51    ReadHistograms(file);
  52 }
  53
  54 AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(const unsigned int harmonic, const bool commonConst, const bool commonHist)
  55    : TNamed(), fHarmonic(harmonic), fNUA(kFALSE), fUseCommonConstants(commonConst), fBookCommonHistograms(commonHist), fCommonHist(0), fQaComponents(0),
  56    fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0),
  57    fPtStatistics(0), fEtaStatistics(0),
  58    fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0)
  59 {
  60    // Constructor defining the harmonic, usage of non uniform acceptance corrections and the AliFlowCommonHist() histograms
  61    // This is the constructor intended for the user
  62    SetNameTitle("MSP","Flow analysis with the Modified Scalar Product method");
  63 }
  64
  65 AliFlowAnalysisWithMSP::~AliFlowAnalysisWithMSP()
  66 {
  67    delete fCommonHist;
  68    delete fQaComponents;
  69    delete fQbComponents;
  70    delete fQaQb;
  71    delete fPtUComponents;
  72    delete fEtaUComponents;
  73    delete fAllStatistics;
  74    delete fPtStatistics;
  75    delete fEtaStatistics;
  76    delete fIntegratedFlow;
  77    delete fDiffFlowPt;
  78    delete fDiffFlowEta;
  79    delete fFlags;
  80 }
  81
  82 void AliFlowAnalysisWithMSP::Init()
  83 {
  84    // Create all output objects. Memory consumption can be reduced by switching off some of the control histograms.
  85    //
  86    delete fCommonHist;     fCommonHist=0;    // Delete existing histograms
  87    delete fQaComponents;   fQaComponents=0;
  88    delete fQbComponents;   fQbComponents=0;
  89    delete fQaQb;           fQaQb=0;
  90    delete fPtUComponents;  fPtUComponents=0;
  91    delete fEtaUComponents; fEtaUComponents=0;
  92    delete fAllStatistics;  fAllStatistics=0;
  93    delete fPtStatistics;   fPtStatistics=0;
  94    delete fEtaStatistics;  fEtaStatistics=0;
  95    delete fIntegratedFlow; fIntegratedFlow=0;
  96    delete fDiffFlowPt;     fDiffFlowPt=0;
  97    delete fDiffFlowEta;    fDiffFlowEta=0;
  98    delete fFlags;          fFlags=0;
  99
 100    // Default binning for histograms
 101    // TODO: allow variable binning
 102
 103    // Defaults
 104    int nBinsPt=100;
 105    double ptMin=0;
 106    double ptMax=10.0;
 107    int nBinsEta=100;
 108    double etaMin=-0.8;
 109    double etaMax=+0.8;
 110
 111    if( fUseCommonConstants || fBookCommonHistograms ) {
 112       AliFlowCommonConstants *c=AliFlowCommonConstants::GetMaster();
 113       nBinsPt=c->GetNbinsPt();
 114       ptMin=c->GetPtMin();
 115       ptMax=c->GetPtMax();
 116       nBinsEta=c->GetNbinsEta();
 117       etaMin=c->GetEtaMin();
 118       etaMax=c->GetEtaMax();
 119    }
 120
 121    if( fBookCommonHistograms ) {           // Use the common constants from the flow package for histogram definitions
 122       std::cerr << "AliFlowAnalysisWithMSP::Init() Creating common histograms" << std::endl;
 123       fCommonHist=new AliFlowCommonHist("AliFlowCommonHist_MSP","AliFlowCommonHist",kTRUE);
 124    }
 125
 126    std::cerr << "AliFlowAnalysisWithMSP::Init() creating MSPHistograms" << std::endl;
 127    // Averages for NUA corrections:
 128    fQaComponents=new AliFlowMSPHistograms(2,"QaComponents",1,0.5,1.5);           // QaX, QaY
 129    fQaComponents->SetXName("all");
 130    fQaComponents->SetVarName("QaX",0);
 131    fQaComponents->SetVarName("QaY",1);
 132    fQbComponents=new AliFlowMSPHistograms(2,"QbComponents",1,0.5,1.5);           // QbX, QbY
 133    fQbComponents->SetXName("all");
 134    fQbComponents->SetVarName("QbX",0);
 135    fQbComponents->SetVarName("QbY",1);
 136    fQaQb=new AliFlowMSPHistograms(1,"QaQb",1,0.5,1.5);                           // QaQb
 137    fQaQb->SetXName("all");
 138    fQaQb->SetVarName("QaQb",0);
 139    fPtUComponents=new AliFlowMSPHistograms(2,"PtUComponents",nBinsPt,ptMin,ptMax);      // ux(pt), uy(pt)
 140    fPtUComponents->SetXName("pt");
 141    fPtUComponents->SetVarName("ux",0);
 142    fPtUComponents->SetVarName("uy",1);
 143    fEtaUComponents=new AliFlowMSPHistograms(2,"EtaUComponents",nBinsEta,etaMin,etaMax);    // ux(eta), uy(eta)
 144    fEtaUComponents->SetXName("eta");
 145    fEtaUComponents->SetVarName("ux",0);
 146    fEtaUComponents->SetVarName("uy",1);
 147
 148    //  Correlation terms
 149    fAllStatistics=new AliFlowMSPHistograms(3,"AllStatistics",1,0.5,1.5);          // terms integrated over pt and eta
 150    fAllStatistics->SetXName("all");
 151    fAllStatistics->SetVarName("uQa",0);
 152    fAllStatistics->SetVarName("uQb",1);
 153    fAllStatistics->SetVarName("QaQb",2);
 154    fPtStatistics=new AliFlowMSPHistograms(3,"PtStatistics",nBinsPt,ptMin,ptMax);         // terms per pt bin
 155    fPtStatistics->SetXName("pt");
 156    fPtStatistics->SetVarName("uQa",0);
 157    fPtStatistics->SetVarName("uQb",1);
 158    fPtStatistics->SetVarName("QaQb",2);
 159    fEtaStatistics=new AliFlowMSPHistograms(3,"EtaStatistics",nBinsEta,etaMin,etaMax);      // terms per eta bin
 160    fEtaStatistics->SetXName("eta");
 161    fEtaStatistics->SetVarName("uQa",0);
 162    fEtaStatistics->SetVarName("uQb",1);
 163    fEtaStatistics->SetVarName("QaQb",2);
 164
 165    fIntegratedFlow=0;   // Created in Finish()
 166    fDiffFlowPt=0;       // Created in Finish
 167    fDiffFlowEta=0;      // Created in Finish
 168
 169    fFlags=new TProfile("Flags","Flags for AliFlowAnalysisWithMSP",10,0.5,10.5,"s");
 170    fFlags->Fill("Harmonic",fHarmonic); // bin 1
 171 }
 172
 173 void AliFlowAnalysisWithMSP::Make(AliFlowEventSimple *event)
 174 {
 175    // Analyze one event. The modified scalar product method estimates flow using the formula:
 176    // BEGIN_LATEX v_2(MSP) = \sqrt( uQ^{a} uQ^{b} / (Q^{a}Q^{b}) ) END_LATEX
 177    // The Q vectors are calculated for the harmonic set previously: fHarmonic.
 178    // Make(event) does not use the new operator, thus avoiding memory leaks in a simple way
 179    // Depending on the compiler about 200 bytes may be pushed/popped on the stack for local variables.
 180    //
 181
 182    if( !event ) return;          // Protect against running without events. Can this happen???
 183
 184    AliFlowVector flowVectors[2];                   // Calculate the two subevent Q vectors:
 185    event->Get2Qsub(flowVectors, fHarmonic);        // No phi, pt or eta weights implemented yet
 186
 187    AliFlowVector &Qa=flowVectors[0];               // Define some mnemonics for the subevent flow vectors:
 188    AliFlowVector &Qb=flowVectors[1];
 189
 190    const double QaW=Qa.GetMult();      // Weight for Qa and combinations of Qa
 191    const double QbW=Qb.GetMult();      // Weight for Qb and combinations of Qb
 192
 193    if( QaW<2 || QbW<2 ) return;  // Require at least 2 particles in each subevent
 194
 195    if(fCommonHist)fCommonHist->FillControlHistograms(event);   // Standard for all flow analysis
 196
 197
 198    const double qaxy[]={Qa.X()/QaW,Qa.Y()/QaW};                // Two variables expected
 199    const double wqaxy[]={QaW,QaW};
 200    fQaComponents->Fill(1, qaxy, wqaxy);                        // only one bin (all pt)
 201
 202    const double qbxy[]={Qb.X()/QbW,Qb.Y()/QbW};                // Two variables expected
 203    const double wqbxy[]={QbW,QbW};
 204    fQbComponents->Fill(1, qbxy, wqbxy);                        // only one bin (all pt)
 205
 206    const double QaQbW=QaW*QbW;
 207    const double weightedQaQb = (Qa*Qb)/QaQbW;                  // Scalar product of subevent Q vectors with weight
 208    fQaQb->Fill(1,&weightedQaQb,&QaQbW);                        // Average of QaQb per event
 209
 210
 211    int iTrack=0;
 212    while( AliFlowTrackSimple *track=event->GetTrack(iTrack++) ) {// Loop over the tracks in the event
 213       // Get the track vector
 214       if(! track->InPOISelection() ) continue;
 215       const double trackWeight=track->Weight();
 216       const double phi=track->Phi();
 217       const double pt=track->Pt();
 218       const double eta=track->Eta();
 219
 220       AliFlowVector u;
 221       u.SetMagPhi(1, fHarmonic*phi, trackWeight);
 222       u.SetMult(1);
 223
 224       // Remove track from subevent a
 225       AliFlowVector mQa(Qa);                 // Initialize with Qa flow vector
 226       if( track->InSubevent(0) ) {
 227          mQa-=u;                             // Should introduce phi weights here
 228       }
 229
 230       // Remove track from subevent b
 231       AliFlowVector mQb(Qb);                 // Initialize with Qb flow vector
 232       if( track->InSubevent(1) ) {
 233          mQb-=u;                             // Should introduce phi weights here
 234       }
 235
 236       const double uQaW = mQa.GetMult();     // Weight is multiplicity of Q vector
 237       const double uQbW = mQb.GetMult();
 238       const double uQa=u*mQa;                // Correlate POI with subevent
 239       const double uQb=u*mQb;
 240
 241       const double uxy[]={u.X(),u.Y()};
 242       const double wxy[]={1.0,1.0};
 243       fPtUComponents->Fill(pt, uxy, wxy);    // vs pt
 244       fEtaUComponents->Fill(eta, uxy, wxy);  // vs eta
 245
 246       const double par[]={uQa/uQaW, uQb/uQbW, weightedQaQb};
 247       const double wgt[]={uQaW, uQbW, QaQbW};
 248       fAllStatistics->Fill(1,   par, wgt  );
 249       fPtStatistics->Fill(pt,   par, wgt );
 250       fEtaStatistics->Fill(eta, par, wgt );
 251    }
 252 }
 253
 254 void AliFlowAnalysisWithMSP::Finish()
 255 {
 256    // Calculate the final result from the stored correlations.
 257    // The NUA corrections are applied if the flag fNUA was set before the call to Finish()
 258    // If the output histograms already exist then they are replaced by the newly calculated result
 259
 260    Print();       // Print a summary of the NUA terms and integrated flow
 261    // TODO: Create result histograms for integrated flow and store the flags to be able to restore the initial state from histograms only
 262
 263    // Create result histograms
 264    fFlags->Fill("NUA",fNUA);
 265    delete fIntegratedFlow; fIntegratedFlow=0;   // First delete existing results (if any)
 266    delete fDiffFlowPt;     fDiffFlowPt=0;
 267    delete fDiffFlowEta;    fDiffFlowEta=0;
 268
 269    fIntegratedFlow=new TH1D("IntegratedFlow","Integrated flow results",10,0.5,10.5);
 270    double vn, vnerror;
 271    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 0);
 272    fIntegratedFlow->SetBinContent(1,vn);
 273    fIntegratedFlow->SetBinError(1,vnerror);
 274    fIntegratedFlow->GetXaxis()->SetBinLabel(1,"POI");
 275    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 1);
 276    fIntegratedFlow->SetBinContent(2,vn);
 277    fIntegratedFlow->SetBinError(2,vnerror);
 278    fIntegratedFlow->GetXaxis()->SetBinLabel(2,"A");
 279    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 2);
 280    fIntegratedFlow->SetBinContent(3,vn);
 281    fIntegratedFlow->SetBinError(3,vnerror);
 282    fIntegratedFlow->GetXaxis()->SetBinLabel(3,"B");
 283
 284    int nbinspt=fPtStatistics->Nbins();
 285    double ptlow=fPtStatistics->XLow();
 286    double pthigh=fPtStatistics->XHigh();
 287    fDiffFlowPt = new TH1D("DiffFlowPt","flow of POI vs pt",nbinspt,ptlow,pthigh);
 288    fDiffFlowPt->SetDirectory(0);          // Do not automatically store in file
 289    fDiffFlowPt->SetStats(kFALSE);
 290    fDiffFlowPt->SetXTitle("p_{t}");
 291    fDiffFlowPt->SetYTitle(Form("v_{%d}",fHarmonic));
 292
 293    for(int i=1; i<=nbinspt; ++i) {
 294       double vn=0;
 295       double evn=0;
 296       if( Calculate(vn, evn, fPtStatistics, fPtUComponents, i) && evn<1 ) {
 297          fDiffFlowPt->SetBinContent(i,vn);
 298          fDiffFlowPt->SetBinError(i,evn);
 299       }
 300    }
 301
 302    int nbinseta=fEtaStatistics->Nbins();
 303    double etalow=fEtaStatistics->XLow();
 304    double etahigh=fEtaStatistics->XHigh();
 305    fDiffFlowEta= new TH1D("DiffFlowEta","flow of POI vs #eta",nbinseta,etalow,etahigh);
 306    fDiffFlowEta->SetDirectory(0);          // Do not automatically store in file
 307    fDiffFlowEta->SetStats(kFALSE);
 308    fDiffFlowEta->SetXTitle("p_{t}");
 309    fDiffFlowEta->SetYTitle(Form("v_{%d}",fHarmonic));
 310
 311    for(int i=1; i<=nbinseta; ++i) {
 312       double vn=0;
 313       double evn=0;
 314       if( Calculate(vn, evn, fEtaStatistics, fEtaUComponents, i) && evn<1 ) {
 315          fDiffFlowEta->SetBinContent(i,vn);
 316          fDiffFlowEta->SetBinError(i,evn);
 317       }
 318    }
 319 }
 320
 321 void AliFlowAnalysisWithMSP::WriteHistograms(TDirectoryFile *file) const
 322 {
 323    //file->Write(file->GetName(), TObject::kSingleKey);      // Make sure the directory itself is written
 324
 325    if(fCommonHist) file->WriteTObject(fCommonHist);
 326
 327    file->WriteTObject(fQaComponents,0,"Overwrite");        // Averages of Qa components per event
 328    file->WriteTObject(fQbComponents,0,"Overwrite");        // Averages of Qb components per event
 329    file->WriteTObject(fQaQb,0,"Overwrite");                // Average of QaQb per event
 330    file->WriteTObject(fPtUComponents,0,"Overwrite");       // u components vs pt
 331    file->WriteTObject(fEtaUComponents,0,"Overwrite");      // u components vs eta
 332    file->WriteTObject(fAllStatistics,0,"Overwrite");       // Integrated uQa, uQb and QaQa
 333    file->WriteTObject(fPtStatistics,0,"Overwrite");        // uQa, uQb and QaQb vs pt
 334    file->WriteTObject(fEtaStatistics,0,"Overwrite");       // uQa, uQb and QaQb vs eta
 335
 336    if( fIntegratedFlow ) file->WriteTObject(fIntegratedFlow,0,"Overwrite");  // Integrated flow for POI and subevents
 337    if( fDiffFlowPt ) file->WriteTObject(fDiffFlowPt,0,"Overwrite");          // Differential flow vs pt if calculated
 338    if( fDiffFlowEta ) file->WriteTObject(fDiffFlowEta,0,"Overwrite");        // Differential flow vs eta if calculated
 339
 340    file->WriteTObject(fFlags,0,"Overwrite");               // fHarmonic, fNUA
 341
 342    file->WriteKeys();   // Make sure it happens now
 343 }
 344
 345 void AliFlowAnalysisWithMSP::WriteCommonResults(TDirectoryFile *file) const
 346 {
 347    // Copy the results to a AliFlowCommonHistResults() class and write to file
 348    // If the results were not calculated again then Finish() is called to generate them
 349    // The global AliFlowCommonConstants() is modified to adapt to the binning used in this analysis
 350
 351    AliFlowCommonConstants *c=AliFlowCommonConstants::GetMaster();
 352    int nBinsPt=fPtStatistics->Nbins();
 353    double ptMin=fPtStatistics->XLow();
 354    double ptMax=fPtStatistics->XHigh();
 355
 356    int nBinsEta=fEtaStatistics->Nbins();
 357    double etaMin=fEtaStatistics->XLow();
 358    double etaMax=fEtaStatistics->XHigh();
 359
 360    c->SetNbinsPt(nBinsPt);
 361    c->SetPtMin(ptMin);
 362    c->SetPtMax(ptMax);
 363    c->SetNbinsEta(nBinsEta);
 364    c->SetEtaMin(etaMin);
 365    c->SetEtaMax(etaMax);
 366
 367    bool oldAddStatus=TH1::AddDirectoryStatus();
 368    TH1::AddDirectory(kFALSE);
 369    AliFlowCommonHistResults *h=new AliFlowCommonHistResults("AliFlowCommonHistResults_MSP","AliFlowCommonHistResults from the MSP method",fHarmonic);
 370
 371    double ivn, ivnerror;
 372    Calculate(ivn, ivnerror, fAllStatistics, fPtUComponents, 0, 0);
 373    h->FillIntegratedFlowPOI(ivn, ivnerror);
 374
 375    for(int bin=1; bin<=nBinsPt; ++bin) {
 376       double vn=0;
 377       double evn=0;
 378       if( Calculate(vn, evn, fPtStatistics, fPtUComponents, bin) && evn>0 ) {
 379          h->FillDifferentialFlowPtPOI(bin, vn, evn);
 380       }
 381    }
 382
 383    for(int bin=1; bin<=nBinsEta; ++bin) {
 384       double vn=0;
 385       double evn=0;
 386       if( Calculate(vn, evn, fEtaStatistics, fEtaUComponents, bin) && evn>0 ) {
 387          h->FillDifferentialFlowEtaPOI(bin, vn, evn);
 388       }
 389    }
 390
 391    file->WriteTObject(h,0,"Overwrite");
 392
 393    TH1::AddDirectory(oldAddStatus);
 394 }
 395
 396
 397
 398 void AliFlowAnalysisWithMSP::Print(const Option_t *opt)const
 399 {
 400    std::cout << setprecision(3);
 401
 402    std::cout << "****************************************************" << std::endl;
 403    std::cout << "    Integrated flow from Modified Scalar Product    " << std::endl;
 404    std::cout << "                                                    " << std::endl;
 405
 406    double vn=0;
 407    double vnerror=0;
 408
 409    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 0);
 410    std::cout << "v" << fHarmonic << " for POI       : " << setw(10) << vn << " +- " << setw(8) << vnerror << std::endl;
 411    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 1);
 412    std::cout << "v" << fHarmonic << " for subevent A: " << setw(10) << vn << " +- " << setw(8) << vnerror << std::endl;
 413    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 2);
 414    std::cout << "v" << fHarmonic << " for subevent B: " << setw(10) << vn << " +- " << setw(8) << vnerror << std::endl;
 415    std::cout << std::endl;
 416
 417    std::cout << "NUA terms: " << (fNUA?"(applied)":"(NOT applied)") << std::endl;
 418
 419    const double ux=fPtUComponents->Average(0);       // Average over all bins
 420    const double eux=TMath::Sqrt(fPtUComponents->Variance(0));
 421    std::cout << "<ux>       " << setw(12) << ux << " +- " << setw(12) << eux << (TMath::Abs(ux)<2*eux?" NOT significant ":" ") << std::endl;
 422    const double ux0eta=fEtaUComponents->Average(fEtaUComponents->FindBin(0.0),0);
 423    const double eux0eta=TMath::Sqrt(fEtaUComponents->Variance(fEtaUComponents->FindBin(0.0),0));
 424    std::cout << "<ux> eta=0 " << setw(12) << ux0eta << " +- " <<  setw(12) << eux0eta << (TMath::Abs(ux0eta)<2*eux0eta?" NOT significant":" ") << std::endl;
 425    const double ux0pt=fPtUComponents->Average(fPtUComponents->FindBin(1.0),0);
 426    const double eux0pt=TMath::Sqrt(fPtUComponents->Variance(fPtUComponents->FindBin(1.0),0));
 427    std::cout << "<ux> pt=1  " << setw(12) << ux0pt << " +- " <<  setw(12) << eux0pt << (TMath::Abs(ux0pt)<2*eux0pt?" NOT significant":" ") << std::endl;
 428
 429    const double uy=fPtUComponents->Average(0);        // Average over all bins
 430    const double euy=TMath::Sqrt(fPtUComponents->Variance(0));
 431    std::cout << "<uy>       " << setw(12) << uy << " +- " << setw(12) << euy << (TMath::Abs(uy)<2*euy?" NOT significant ":" ") << std::endl;;
 432    const double uy0eta=fEtaUComponents->Average(fEtaUComponents->FindBin(0.0),1);
 433    const double euy0eta=TMath::Sqrt(fEtaUComponents->Variance(fEtaUComponents->FindBin(0.0),1));
 434    std::cout << "<uy> eta=0 " << setw(12) << uy0eta << " +- " << setw(12) << euy0eta << (TMath::Abs(uy0eta)<2*euy0eta?" NOT significant ":" ") << std::endl;
 435    const double uy0pt=fPtUComponents->Average(fPtUComponents->FindBin(1.0),1);
 436    const double euy0pt=TMath::Sqrt(fPtUComponents->Variance(fPtUComponents->FindBin(1.0),1));
 437    std::cout << "<uy> pt=1  " << setw(12) << uy0pt << " +- " << setw(12) << euy0pt << (TMath::Abs(uy0pt)<2*euy0pt?" NOT significant ":" ") << std::endl;
 438
 439    std::cout << "<QaX>      " << setw(12) << fQaComponents->Average(0) << " +- " << setw(12) << TMath::Sqrt(fQaComponents->Variance(0)) << std::endl;
 440    std::cout << "<QaY>      " << setw(12) << fQaComponents->Average(1) << " +- " << setw(12) << TMath::Sqrt(fQaComponents->Variance(1)) << std::endl;
 441    std::cout << "<QbX>      " << setw(12) << fQbComponents->Average(0) << " +- " << setw(12) << TMath::Sqrt(fQbComponents->Variance(0)) << std::endl;
 442    std::cout << "<QbY>      " << setw(12) << fQbComponents->Average(1) << " +- " << setw(12) << TMath::Sqrt(fQbComponents->Variance(1)) << std::endl;
 443    std::cout << "<QaQb>     " << setw(12) << fQaQb->Average(0) << " +- " << setw(12) << TMath::Sqrt(fQbComponents->Variance(0)) << std::endl;
 444    std::cout << std::endl;
 445
 446    std::cout << "Covariance matrix: " << std::endl;
 447    std::cout << "     " << setw(12) << "uQa"                           << setw(12) << "uQb"                           << setw(12) << "QaQb" << std::endl;
 448    std::cout << "uQa  " << setw(12) << fAllStatistics->Covariance(0,0) << std::endl;
 449    std::cout << "uQb  " << setw(12) << fAllStatistics->Covariance(1,0) << setw(12) << fAllStatistics->Covariance(1,1) << std::endl;
 450    std::cout << "QaQb " << setw(12) << fAllStatistics->Covariance(2,0) << setw(12) << fAllStatistics->Covariance(2,1) << setw(12) << fQaQb->Variance(0) << std::endl;
 451    std::cout << "****************************************************" << std::endl;
 452    std::cout << std::endl;
 453 }
 454
 455
 456 // private functions --------------------------------------------------------------------------------------
 457 bool AliFlowAnalysisWithMSP::Calculate(double &vn, double &vnerror, const AliFlowMSPHistograms *hist, const AliFlowMSPHistograms *components, const int bin, const int poi) const
 458 {
 459    // Get all averages and correlations need for the flow calculation
 460    double uQa=hist->Average(bin,0);                // <<uQa>>
 461    double VuQa=hist->Variance(bin,0);              // Var(<<uQa>>)
 462    double uQb=hist->Average(bin,1);                // <<uQb>>
 463    double VuQb=hist->Variance(bin,1);              // Var(<<uQb>>)
 464    double QaQb=fQaQb->Average(1,0);                // <QaQb> Should not be taken from hist(bin) buit from fQaQa because there QaQb is entered multiple times: <<QaQb>>!!
 465    double VQaQb=fQaQb->Variance(1,0);              // V(<QaQb>)
 466    double CuQauQb=hist->Covariance(bin,0,1);       // Cov(<<uQa>>,<<uQb>>)
 467    double CuQaQaQb=hist->Covariance(bin,0,2);      // Cov(<<uQa>>,<QaQb>)
 468    double CuQbQaQb=hist->Covariance(bin,1,2);      // Cov(<<uQb>>,<QaQb>)
 469
 470    if( fNUA && components ) {
 471       // Apply NUA correction to QaQb.
 472       double QaX=fQaComponents->Average(0);
 473       double QaY=fQaComponents->Average(1);
 474       double QbX=fQbComponents->Average(0);
 475       double QbY=fQbComponents->Average(1);
 476
 477       QaQb=QaQb-QaX*QbX-QaY*QbY;
 478
 479       // Apply NUA to uQa and uQb (per bin)
 480       double uX=components->Average(bin,0);  // bin 0 is integrated over all bins
 481       double uY=components->Average(bin,1);
 482
 483       uQa = uQa - uX*QaX - uY*QaY;
 484       uQb = uQb - uX*QbX - uY*QbY;
 485       // Error calculation not fully modified: only above terms, the spread in <<u>> and <Qa> and <Qb> are not used
 486       // therefore this should only be applied if significant!
 487       // TODO: Check if not fully NUA correcting the error calculation is justified (but this is compatible with the original SP method)!
 488    }
 489
 490    // Some sanity checks:
 491    if( uQa*uQb*QaQb <= 0 ) {                                    // Catch imaginary results
 492       vn=-99;
 493       vnerror=-99;
 494       return false;
 495    }
 496
 497    // Decent results, calculate, print and store
 498    // TODO: The three cases are cyclic permutations of u, Qa, Qb so there should be a simpler way to do this.
 499    // However the difficulty is that we deal here with uQa, uQb and QaQb
 500    switch (poi) {
 501    case 1:                                        // Subevent A reference flow
 502       {
 503          if( TMath::Abs(uQb) < 1e-30*TMath::Abs(uQa*QaQb) ) {  // Protect against infinity
 504             vn=0;
 505             vnerror=-1;
 506             return false;
 507          }
 508          double vnB = TMath::Sqrt( uQa*QaQb / uQb ); // vn
 509          double VvnB = QaQb*VuQa/(4*uQa*uQb)         // Variance of vn
 510             + uQa*VQaQb/(4*QaQb*uQb)
 511             + uQa*QaQb*VuQb/(4*TMath::Power(uQb,3))
 512             + CuQaQaQb/(2*uQb)
 513             - QaQb*CuQauQb/(2*TMath::Power(uQb,2))
 514             - uQa*CuQbQaQb/(2*TMath::Power(uQb,2));
 515          vn=vnB;
 516          if( VvnB<0 ) {
 517             vnerror=VvnB;
 518             return false;
 519          }
 520          vnerror=TMath::Sqrt(VvnB);
 521       }
 522       break;
 523    case 2:                                        // Subevent B reference flow
 524       {
 525          if( TMath::Abs(uQa) < 1e-30*TMath::Abs(uQb*QaQb) ) {  // Protect against infinity
 526             vn=0;
 527             vnerror=-1;
 528             return false;
 529          }
 530          double vnA = TMath::Sqrt( uQb*QaQb / uQa );        // vn
 531          double VvnA = uQb*VQaQb/(4*QaQb*uQa)               // Variance of vn
 532             + QaQb*VuQb/(4*uQb*uQa)
 533             + QaQb*uQb*VuQa/(4*TMath::Power(uQa,3))
 534             + CuQbQaQb/(2*uQa)
 535             - uQb*CuQaQaQb/(2*TMath::Power(uQa,2))
 536             - uQa*CuQauQb/(2*TMath::Power(uQa,2));
 537          vn=vnA;
 538          if( VvnA<0 ) {
 539             vnerror=VvnA;
 540             return false;
 541          }
 542          vnerror=TMath::Sqrt(VvnA);
 543       }
 544       break;
 545    default:                                       // POI flow
 546       {
 547          if( TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQa*uQb) ) {  // Catch infinity
 548             vn=0;
 549             vnerror=-1;
 550             return false;
 551          }
 552          double vnP = TMath::Sqrt( uQa*uQb / QaQb );        // vn
 553          if(   TMath::Abs(uQa*QaQb) < 1e-30*TMath::Abs(uQb*VuQa)
 554             || TMath::Abs(uQb*QaQb) < 1e-30*TMath::Abs(uQa*VuQb)
 555             || TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQa*uQb*VQaQb)
 556             || TMath::Abs(QaQb) < 1e-30*TMath::Abs(CuQauQb)
 557             || TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQb*CuQaQaQb)
 558             || TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQa*CuQbQaQb)
 559             ){
 560                vnerror=-98;
 561                return false;
 562          }
 563          double VvnP = uQb*VuQa/(4*uQa*QaQb)                // Variance of vn
 564             + uQa*VuQb/(4*uQb*QaQb)
 565             + uQa*uQb*VQaQb/(4*TMath::Power(QaQb,3))
 566             + CuQauQb/(2*QaQb)
 567             - uQb*CuQaQaQb/(2*TMath::Power(QaQb,2))
 568             - uQa*CuQbQaQb/(2*TMath::Power(QaQb,2));
 569          vn=vnP;
 570          if( VvnP<0 ) {
 571             vnerror=VvnP;
 572             return false;
 573          }
 574          vnerror=TMath::Sqrt(VvnP);
 575       }
 576    }  // Switch between POI and subevents
 577
 578    return (vnerror>=0);
 579 }
 580
 581 void AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *file)
 582 {
 583    delete fCommonHist;     fCommonHist=0;    // Delete existing histograms
 584    delete fQaComponents;   fQaComponents=0;
 585    delete fQbComponents;   fQbComponents=0;
 586    delete fQaQb;           fQaQb=0;
 587    delete fPtUComponents;  fPtUComponents=0;
 588    delete fEtaUComponents; fEtaUComponents=0;
 589    delete fAllStatistics;  fAllStatistics=0;
 590    delete fPtStatistics;   fPtStatistics=0;
 591    delete fEtaStatistics;  fEtaStatistics=0;
 592    delete fIntegratedFlow; fIntegratedFlow=0;
 593    delete fDiffFlowPt;     fDiffFlowPt=0;
 594    delete fDiffFlowEta;    fDiffFlowEta=0;
 595    delete fFlags;          fFlags=0;
 596
 597    file->GetObject("QaComponents",fQaComponents);
 598    file->GetObject("QbComponents",fQbComponents);
 599    file->GetObject("QaQb",fQaQb);
 600    file->GetObject("PtUComponents",fPtUComponents);
 601    file->GetObject("EtaUComponents",fEtaUComponents);
 602    file->GetObject("AllStatistics",fAllStatistics);
 603    file->GetObject("PtStatistics",fPtStatistics);
 604    file->GetObject("EtaStatistics",fEtaStatistics);
 605    if( !fQaComponents || !fQbComponents || !fQaQb || !fPtUComponents || !fEtaUComponents || !fAllStatistics || !fPtStatistics || !fEtaStatistics ) {
 606       std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms() : One or more histograms were not read correctly from " << file->GetPath() << std::endl;
 607    }
 608
 609    // Optional histograms, may return a zero pointer
 610    file->GetObject("AliFlowCommonHist_MSP",fCommonHist);              // The AliFlowCommonHist is optional
 611    file->GetObject("IntegratedFlow",fIntegratedFlow);
 612    file->GetObject("DiffFlowPt",fDiffFlowPt);
 613    file->GetObject("DiffFlowEta",fDiffFlowEta);
 614
 615    file->GetObject("Flags",fFlags);
 616
 617    if( !fFlags ){
 618       std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms() : Flags histogram not found, using defaults" << std::endl;
 619       fHarmonic=2;
 620       fNUA=false;
 621    }else{
 622       fHarmonic=fFlags->GetBinContent(1);
 623       double harmonicSpread=fFlags->GetBinError(1);
 624       if( harmonicSpread!=0 ) {
 625          std::cerr << "These histograms seem to be merged from analysis with different Harmonics. Results are invalid!" << std::endl;
 626          delete fIntegratedFlow; fIntegratedFlow=0;
 627          delete fDiffFlowPt;     fDiffFlowPt=0;
 628          delete fDiffFlowEta;    fDiffFlowEta=0;
 629       }
 630       fNUA=fFlags->GetBinContent(2);   // Mixing NUA does not matter since it needs to be recalculated anyway
 631       double nuaSpread=fFlags->GetBinError(2);
 632       if( nuaSpread!=0 ) {
 633          std::cerr << "These histograms seem to be merged from analysis with different NUA corrections. Results removed" << std::endl;
 634          delete fIntegratedFlow; fIntegratedFlow=0;
 635          delete fDiffFlowPt;     fDiffFlowPt=0;
 636          delete fDiffFlowEta;    fDiffFlowEta=0;
 637       }
 638    }
 639    fBookCommonHistograms=(fCommonHist!=0);
 640 }