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(),
  38    fHarmonic(2), fNUA(kFALSE), fUseCommonConstants(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0),
  39    fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0),
  40    fPtStatistics(0), fEtaStatistics(0),
  41    fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0), fHistList(0)
  42 {
  43    // Default constructor. Intended for root IO purposes only
  44 }
  45
  46 AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(TDirectory *file)
  47    : TNamed(),
  48    fHarmonic(2), fNUA(kFALSE), fUseCommonConstants(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0),
  49    fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0),
  50    fPtStatistics(0), fEtaStatistics(0),
  51    fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0), fHistList(0)
  52 {
  53    // Constructor reads the internal state from a root file.
  54    // No check for consistency is done. If flags or histograms are not found they are left at default (0).
  55    ReadHistograms(file);
  56 }
  57
  58 AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(TList *list)
  59    : TNamed(),
  60    fHarmonic(2), fNUA(kFALSE), fUseCommonConstants(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0),
  61    fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0),
  62    fPtStatistics(0), fEtaStatistics(0),
  63    fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0), fHistList(0)
  64 {
  65    // Constructor reads the internal state from a root file.
  66    // No check for consistency is done. If flags or histograms are not found they are left at default (0).
  67    ReadHistograms(list);
  68 }
  69
  70
  71 AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(const unsigned int harmonic, const bool commonConst, const bool commonHist)
  72    : TNamed(),
  73    fHarmonic(harmonic), fNUA(kFALSE), fUseCommonConstants(commonConst), fBookCommonHistograms(commonHist), fCommonHist(0), fQaComponents(0),
  74    fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0),
  75    fPtStatistics(0), fEtaStatistics(0),
  76    fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0), fHistList(0)
  77 {
  78    // Constructor defining the harmonic, usage of non uniform acceptance corrections and the AliFlowCommonHist() histograms
  79    // Calls Init() to set-up the histograms. This is equivalent to:
  80    // AliFlowAnalysisWithMSP *ana= new AliFlowAnalysisWithMSP);
  81    // ana->SetHarmonic(harmonic);
  82    // ana->UseCommonConstants(commonConst);
  83    // ana->EnableCommonHistograms(commonHist);
  84    // ana->Init();
  85    // This is the constructor intended for the user
  86    SetNameTitle("MSP","Flow analysis with the Modified Scalar Product method");
  87    Init();
  88 }
  89
  90 AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(const AliFlowAnalysisWithMSP &x)
  91    : TNamed(),
  92    fHarmonic(x.fHarmonic), fNUA(x.fNUA), fUseCommonConstants(x.fUseCommonConstants), fBookCommonHistograms(x.fBookCommonHistograms), fCommonHist(0),
  93    fQaComponents(0), fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0),
  94    fPtStatistics(0), fEtaStatistics(0),
  95    fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0), fHistList(0)
  96 {
  97    // Copy constructor
  98    SetNameTitle("MSP","Flow analysis with the Modified Scalar Product method");
  99    if( x.fQaComponents )   fQaComponents=(AliFlowMSPHistograms *)(x.fQaComponents)->Clone();
 100    if( x.fQbComponents )   fQbComponents=(AliFlowMSPHistograms *)(x.fQbComponents)->Clone();
 101    if( x.fQaQb )           fQaQb=(AliFlowMSPHistograms *)(x.fQaQb)->Clone();
 102    if( x.fPtUComponents)     fPtUComponents=(AliFlowMSPHistograms *)(x.fPtUComponents)->Clone();
 103    if( x.fEtaUComponents )   fEtaUComponents=(AliFlowMSPHistograms *)(x.fEtaUComponents)->Clone();
 104    if( x.fAllStatistics )    fAllStatistics=(AliFlowMSPHistograms *)(x.fAllStatistics)->Clone();
 105    if( x.fPtStatistics )     fPtStatistics=(AliFlowMSPHistograms *)(x.fPtStatistics)->Clone();
 106    if( x.fEtaStatistics )    fEtaStatistics=(AliFlowMSPHistograms *)(x.fEtaStatistics)->Clone();
 107    if( x.fIntegratedFlow )   fIntegratedFlow=(TH1D *)(x.fIntegratedFlow)->Clone();
 108    if( x.fDiffFlowPt )       fDiffFlowPt=(TH1D *)(x.fDiffFlowPt)->Clone();
 109    if( x.fDiffFlowEta )      fDiffFlowEta=(TH1D *)(x.fDiffFlowEta)->Clone();
 110    if( x.fFlags )            fFlags=(TProfile *)(x.fFlags)->Clone();
 111    if( x.fCommonHist )       fCommonHist=new AliFlowCommonHist(*(x.fCommonHist));
 112    // The histogram list fHistList is not cloned because it is regenerated when requested
 113 }
 114
 115 AliFlowAnalysisWithMSP::~AliFlowAnalysisWithMSP()
 116 {
 117    // Destructor. All internal objects are owned by this class and deleted here.
 118    delete fCommonHist;
 119    delete fQaComponents;
 120    delete fQbComponents;
 121    delete fQaQb;
 122    delete fPtUComponents;
 123    delete fEtaUComponents;
 124    delete fAllStatistics;
 125    delete fPtStatistics;
 126    delete fEtaStatistics;
 127    delete fIntegratedFlow;
 128    delete fDiffFlowPt;
 129    delete fDiffFlowEta;
 130    delete fFlags;
 131    if(fHistList) {
 132       fHistList->SetOwner(kFALSE);   // Histograms were already deleted manually
 133       delete fHistList;              // Delete the TList itself
 134    }
 135 }
 136
 137 void AliFlowAnalysisWithMSP::Init()
 138 {
 139    // Create all output objects. Memory consumption can be reduced by switching off some of the control histograms.
 140    // pt and eta binning are set to the values defined in the static class AliFlowCommonConstants if enabled. Otherwise defaults are set.
 141    delete fCommonHist;     fCommonHist=0;    // Delete existing histograms
 142    delete fQaComponents;   fQaComponents=0;
 143    delete fQbComponents;   fQbComponents=0;
 144    delete fQaQb;           fQaQb=0;
 145    delete fPtUComponents;  fPtUComponents=0;
 146    delete fEtaUComponents; fEtaUComponents=0;
 147    delete fAllStatistics;  fAllStatistics=0;
 148    delete fPtStatistics;   fPtStatistics=0;
 149    delete fEtaStatistics;  fEtaStatistics=0;
 150    delete fIntegratedFlow; fIntegratedFlow=0;
 151    delete fDiffFlowPt;     fDiffFlowPt=0;
 152    delete fDiffFlowEta;    fDiffFlowEta=0;
 153    delete fFlags;          fFlags=0;
 154    if(fHistList) {
 155       fHistList->SetOwner(kFALSE);   // Histograms were already deleted manually
 156       delete fHistList;              // Delete the TList itself
 157       fHistList=0;                   // Clear pointer which is invalid afcter delete
 158    }
 159
 160    // Default binning for histograms
 161    // TODO: allow variable binning in a simpler way than by AliFlowCommonConstants() only
 162
 163    // Defaults
 164    int nBinsPt=100;
 165    double ptMin=0;
 166    double ptMax=10.0;
 167    int nBinsEta=100;
 168    double etaMin=-0.8;
 169    double etaMax=+0.8;
 170
 171    if( fUseCommonConstants || fBookCommonHistograms ) {
 172       AliFlowCommonConstants *c=AliFlowCommonConstants::GetMaster();
 173       nBinsPt=c->GetNbinsPt();
 174       ptMin=c->GetPtMin();
 175       ptMax=c->GetPtMax();
 176       nBinsEta=c->GetNbinsEta();
 177       etaMin=c->GetEtaMin();
 178       etaMax=c->GetEtaMax();
 179    }
 180
 181    if( fBookCommonHistograms ) {           // Use the common constants from the flow package for histogram definitions
 182       std::cerr << "AliFlowAnalysisWithMSP::Init() Creating common histograms" << std::endl;
 183       fCommonHist=new AliFlowCommonHist("AliFlowCommonHist_MSP","AliFlowCommonHist",kTRUE);
 184    }
 185
 186    std::cerr << "AliFlowAnalysisWithMSP::Init() creating MSPHistograms" << std::endl;
 187    // Averages for NUA corrections:
 188    fQaComponents=new AliFlowMSPHistograms(2,"QaComponents",1,0.5,1.5);           // QaX, QaY
 189    fQaComponents->SetXName("all");
 190    fQaComponents->SetVarName("QaX",0);
 191    fQaComponents->SetVarName("QaY",1);
 192    fQbComponents=new AliFlowMSPHistograms(2,"QbComponents",1,0.5,1.5);           // QbX, QbY
 193    fQbComponents->SetXName("all");
 194    fQbComponents->SetVarName("QbX",0);
 195    fQbComponents->SetVarName("QbY",1);
 196    fQaQb=new AliFlowMSPHistograms(1,"QaQb",1,0.5,1.5);                           // QaQb
 197    fQaQb->SetXName("all");
 198    fQaQb->SetVarName("QaQb",0);
 199    fPtUComponents=new AliFlowMSPHistograms(2,"PtUComponents",nBinsPt,ptMin,ptMax);      // ux(pt), uy(pt)
 200    fPtUComponents->SetXName("pt");
 201    fPtUComponents->SetVarName("ux",0);
 202    fPtUComponents->SetVarName("uy",1);
 203    fEtaUComponents=new AliFlowMSPHistograms(2,"EtaUComponents",nBinsEta,etaMin,etaMax);    // ux(eta), uy(eta)
 204    fEtaUComponents->SetXName("eta");
 205    fEtaUComponents->SetVarName("ux",0);
 206    fEtaUComponents->SetVarName("uy",1);
 207
 208    //  Correlation terms
 209    fAllStatistics=new AliFlowMSPHistograms(3,"AllStatistics",1,0.5,1.5);          // terms integrated over pt and eta
 210    fAllStatistics->SetXName("all");
 211    fAllStatistics->SetVarName("uQa",0);
 212    fAllStatistics->SetVarName("uQb",1);
 213    fAllStatistics->SetVarName("QaQb",2);
 214    fPtStatistics=new AliFlowMSPHistograms(3,"PtStatistics",nBinsPt,ptMin,ptMax);         // terms per pt bin
 215    fPtStatistics->SetXName("pt");
 216    fPtStatistics->SetVarName("uQa",0);
 217    fPtStatistics->SetVarName("uQb",1);
 218    fPtStatistics->SetVarName("QaQb",2);
 219    fEtaStatistics=new AliFlowMSPHistograms(3,"EtaStatistics",nBinsEta,etaMin,etaMax);      // terms per eta bin
 220    fEtaStatistics->SetXName("eta");
 221    fEtaStatistics->SetVarName("uQa",0);
 222    fEtaStatistics->SetVarName("uQb",1);
 223    fEtaStatistics->SetVarName("QaQb",2);
 224
 225    fIntegratedFlow=0;   // Created in Finish()
 226    fDiffFlowPt=0;       // Created in Finish
 227    fDiffFlowEta=0;      // Created in Finish
 228
 229    fFlags=new TProfile("Flags","Flags for AliFlowAnalysisWithMSP",10,0.5,10.5,"s");
 230    fFlags->Fill("Harmonic",fHarmonic); // bin 1
 231 }
 232
 233 void AliFlowAnalysisWithMSP::Make(AliFlowEventSimple *event)
 234 {
 235    // Analyze one event. The modified scalar product method estimates flow using the formula:
 236    // BEGIN_LATEX v_2(MSP) = \sqrt( uQ^{a} uQ^{b} / (Q^{a}Q^{b}) ) END_LATEX
 237    // The Q vectors are calculated for the harmonic set previously: fHarmonic.
 238    // Make(event) does not use the new operator, thus avoiding memory leaks in a simple way
 239    // Depending on the compiler about 200 bytes may be pushed/popped on the stack for local variables.
 240    //
 241
 242    if( !event ) return;          // Protect against running without events. Can this happen???
 243
 244    AliFlowVector flowVectors[2];                   // Calculate the two subevent Q vectors:
 245    event->Get2Qsub(flowVectors, fHarmonic);        // TODO: Check how do the phi, pt, eta weights work in the flow event?
 246
 247    AliFlowVector &Qa=flowVectors[0];               // Define some mnemonics for the subevent flow vectors:
 248    AliFlowVector &Qb=flowVectors[1];
 249
 250    const double QaW=Qa.GetMult();      // Weight for Qa and combinations of Qa
 251    const double QbW=Qb.GetMult();      // Weight for Qb and combinations of Qb
 252
 253    if( QaW<2 || QbW<2 ) return;  // Require at least 2 particles in each subevent
 254
 255    if(fCommonHist)fCommonHist->FillControlHistograms(event);   // Standard for all flow analysis
 256
 257    // Fill NUA correction histograms for Qa, Qb and QaQb per event:
 258    const double qaxy[]={Qa.X()/QaW,Qa.Y()/QaW};                // Two variables: Qa components
 259    const double wqaxy[]={QaW,QaW};
 260    fQaComponents->Fill(1, qaxy, wqaxy);                        // only one bin (all pt)
 261
 262    const double qbxy[]={Qb.X()/QbW,Qb.Y()/QbW};                // Two variables: Qb components
 263    const double wqbxy[]={QbW,QbW};
 264    fQbComponents->Fill(1, qbxy, wqbxy);                        // only one bin (all pt)
 265
 266    const double QaQbW[]={QaW*QbW};
 267    const double weightedQaQb[] = {(Qa*Qb)/QaQbW[0]};                  // Scalar product of subevent Q vectors with weight, only 1 variable
 268    fQaQb->Fill(1,weightedQaQb,QaQbW);                        // Average of QaQb per event
 269
 270    int iTrack=0;
 271    while( AliFlowTrackSimple *track=event->GetTrack(iTrack++) ) {// Loop over the tracks in the event
 272       if(! track->InPOISelection() ) continue;  // Ignore if not a POI
 273       const double trackWeight=track->Weight(); // Get the track vector
 274       const double phi=track->Phi();
 275       const double pt=track->Pt();
 276       const double eta=track->Eta();
 277
 278       AliFlowVector u;
 279       u.SetMagPhi(trackWeight, fHarmonic*phi, trackWeight);    // Length = track weight = multiplicity for a single track
 280
 281       // Remove track from subevent a
 282       AliFlowVector mQa(Qa);                 // Initialize with Qa flow vector of this event
 283       if( track->InSubevent(0) ) {           // Correct for autocorrelation with POI for this track
 284          mQa-=u;
 285       }
 286
 287       // Remove track from subevent b
 288       AliFlowVector mQb(Qb);                 // Initialize with Qb flow vector of this event
 289       if( track->InSubevent(1) ) {           // Correct for autocorrelation with POI for this track
 290          mQb-=u;
 291       }
 292
 293       const double uQaW = mQa.GetMult()*trackWeight;     // Weight is multiplicity of Q vector times weight of POI
 294       const double uQbW = mQb.GetMult()*trackWeight;
 295       const double uQa=u*mQa;                            // Correlate POI with subevent
 296       const double uQb=u*mQb;
 297
 298       const double uxy[]={u.X(),u.Y()};
 299       const double wxy[]={1.0,1.0};
 300       fPtUComponents->Fill(pt, uxy, wxy);    // NUA for POI vs pt
 301       fEtaUComponents->Fill(eta, uxy, wxy);  // NUA for POI vs eta
 302
 303       const double par[]={uQa/uQaW, uQb/uQbW, weightedQaQb[0]};   // Three variables to correlate: uQa, uQb and QaQb
 304       const double wgt[]={uQaW, uQbW, QaQbW[0]};
 305       fAllStatistics->Fill(1,   par, wgt  );    // only 1 bin, integrated over pt and eta
 306       fPtStatistics->Fill(pt,   par, wgt );     // pt differential correlations
 307       fEtaStatistics->Fill(eta, par, wgt );     // eta differential correlations
 308    }
 309 }
 310
 311 void AliFlowAnalysisWithMSP::Finish()
 312 {
 313    // Calculate the final result from the stored correlations.
 314    // The NUA corrections are applied if the flag fNUA was set before the call to Finish()
 315    // If the output histograms already exist then they are replaced by the newly calculated result
 316
 317    Print();       // Print a summary of the NUA terms and integrated flow
 318
 319    // Create result histograms
 320    fFlags->Fill("NUA",fNUA);
 321    delete fIntegratedFlow; fIntegratedFlow=0;   // First delete existing results (if any)
 322    delete fDiffFlowPt;     fDiffFlowPt=0;
 323    delete fDiffFlowEta;    fDiffFlowEta=0;
 324
 325    bool oldAddStatus=TH1::AddDirectoryStatus();       // Do not store the next histograms automatically
 326    TH1::AddDirectory(kFALSE);                         // We need full control over the writing of the hostograms
 327
 328    fIntegratedFlow=new TH1D("IntegratedFlow","Integrated flow results",10,0.5,10.5);
 329    fIntegratedFlow->SetDirectory(0);
 330    double vn, vnerror;
 331    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 0);
 332    fIntegratedFlow->SetBinContent(1,vn);
 333    fIntegratedFlow->SetBinError(1,vnerror);
 334    fIntegratedFlow->GetXaxis()->SetBinLabel(1,"POI");
 335    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 1);
 336    fIntegratedFlow->SetBinContent(2,vn);
 337    fIntegratedFlow->SetBinError(2,vnerror);
 338    fIntegratedFlow->GetXaxis()->SetBinLabel(2,"A");
 339    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 2);
 340    fIntegratedFlow->SetBinContent(3,vn);
 341    fIntegratedFlow->SetBinError(3,vnerror);
 342    fIntegratedFlow->GetXaxis()->SetBinLabel(3,"B");
 343
 344    int nbinspt=fPtStatistics->Nbins();
 345    double ptlow=fPtStatistics->XLow();
 346    double pthigh=fPtStatistics->XHigh();
 347    fDiffFlowPt = new TH1D("DiffFlowPt","flow of POI vs pt",nbinspt,ptlow,pthigh);
 348    fDiffFlowPt->SetDirectory(0);          // Do not automatically store in file
 349    fDiffFlowPt->SetStats(kFALSE);
 350    fDiffFlowPt->SetXTitle("p_{t}");
 351    fDiffFlowPt->SetYTitle(Form("v_{%d}",fHarmonic));
 352
 353    for(int i=1; i<=nbinspt; ++i) {
 354       double vnpt=0;
 355       double evnpt=0;
 356       if( Calculate(vnpt, evnpt, fPtStatistics, fPtUComponents, i) && evnpt<1 ) {
 357          fDiffFlowPt->SetBinContent(i,vnpt);
 358          fDiffFlowPt->SetBinError(i,evnpt);
 359       }
 360    }
 361
 362    int nbinseta=fEtaStatistics->Nbins();
 363    double etalow=fEtaStatistics->XLow();
 364    double etahigh=fEtaStatistics->XHigh();
 365    fDiffFlowEta= new TH1D("DiffFlowEta","flow of POI vs #eta",nbinseta,etalow,etahigh);
 366    fDiffFlowEta->SetDirectory(0);          // Do not automatically store in file
 367    fDiffFlowEta->SetStats(kFALSE);
 368    fDiffFlowEta->SetXTitle("#eta");
 369    fDiffFlowEta->SetYTitle(Form("v_{%d}",fHarmonic));
 370
 371    for(int i=1; i<=nbinseta; ++i) {
 372       double vneta=0;
 373       double evneta=0;
 374       if( Calculate(vneta, evneta, fEtaStatistics, fEtaUComponents, i) && evneta<1 ) {
 375          fDiffFlowEta->SetBinContent(i,vneta);
 376          fDiffFlowEta->SetBinError(i,evneta);
 377       }
 378    }
 379
 380    TH1::AddDirectory(oldAddStatus);
 381 }
 382
 383 void AliFlowAnalysisWithMSP::WriteHistograms(TDirectoryFile *file) const
 384 {
 385    // Write the internal status to file. If the AliFlowCommonHistograms were enabled they are written also.
 386    // Results are written only if they were calculated by Finish().
 387    // From these histograms the internal state of *this can be reconstructed with ReadHistograms()
 388
 389    //file->Write(file->GetName(), TObject::kSingleKey);     // Make sure the directory itself is written
 390
 391    if(fCommonHist) file->WriteTObject(fCommonHist);
 392
 393    TList *t=new TList();                                    // This object is written as a marker for redoFinish()
 394    t->SetName("cobjMSP");
 395    file->WriteTObject(t,0,"Overwrite");
 396    file->WriteTObject(fQaComponents,0,"Overwrite");         // Averages of Qa components per event
 397    file->WriteTObject(fQbComponents,0,"Overwrite");         // Averages of Qb components per event
 398    file->WriteTObject(fQaQb,0,"Overwrite");                 // Average of QaQb per event
 399    file->WriteTObject(fPtUComponents,0,"Overwrite");        // u components vs pt
 400    file->WriteTObject(fEtaUComponents,0,"Overwrite");       // u components vs eta
 401    file->WriteTObject(fAllStatistics,0,"Overwrite");        // Integrated uQa, uQb and QaQa
 402    file->WriteTObject(fPtStatistics,0,"Overwrite");         // uQa, uQb and QaQb vs pt
 403    file->WriteTObject(fEtaStatistics,0,"Overwrite");        // uQa, uQb and QaQb vs eta
 404
 405    if( fIntegratedFlow ) file->WriteTObject(fIntegratedFlow,0,"Overwrite");  // Integrated flow for POI and subevents
 406    if( fDiffFlowPt ) file->WriteTObject(fDiffFlowPt,0,"Overwrite");          // Differential flow vs pt if calculated
 407    if( fDiffFlowEta ) file->WriteTObject(fDiffFlowEta,0,"Overwrite");        // Differential flow vs eta if calculated
 408
 409    file->WriteTObject(fFlags,0,"Overwrite");               // fHarmonic, fNUA
 410
 411    file->WriteKeys();   // Make sure it happens now
 412 }
 413
 414 void AliFlowAnalysisWithMSP::WriteCommonResults(TDirectoryFile *file) const
 415 {
 416    // Export the results to a AliFlowCommonHistResults() class and write to file
 417    // If the results were not calculated then Finish() is called to generate them
 418
 419    int nBinsPt=fPtStatistics->Nbins();       // Get the actually used binning from the Statistics histograms
 420    double ptMin=fPtStatistics->XLow();
 421    double ptMax=fPtStatistics->XHigh();
 422
 423    int nBinsEta=fEtaStatistics->Nbins();
 424    double etaMin=fEtaStatistics->XLow();
 425    double etaMax=fEtaStatistics->XHigh();
 426
 427    AliFlowCommonConstants *c=AliFlowCommonConstants::GetMaster(); // Get the static common constants object
 428    // Save the old AliFlowCommonConstants status
 429    int oldNbinsPt=c->GetNbinsPt();
 430    double oldPtMin=c->GetPtMin();
 431    double oldPtMax=c->GetPtMax();
 432    int oldNbinsEta=c->GetNbinsEta();
 433    double oldEtaMin=c->GetEtaMin();
 434    double oldEtaMax=c->GetEtaMax();
 435    // Modify AliFlowCommonConstants to make sure that AliFlowCommonResults is generated with the correct binning
 436    c->SetNbinsPt(nBinsPt);
 437    c->SetPtMin(ptMin);
 438    c->SetPtMax(ptMax);
 439    c->SetNbinsEta(nBinsEta);
 440    c->SetEtaMin(etaMin);
 441    c->SetEtaMax(etaMax);
 442
 443    bool oldAddStatus=TH1::AddDirectoryStatus();       // Do not store the next histograms automatically
 444    TH1::AddDirectory(kFALSE);                         // We need full control over the writing of the hostograms
 445    AliFlowCommonHistResults *h=new AliFlowCommonHistResults("AliFlowCommonHistResults_MSP","AliFlowCommonHistResults from the MSP method",fHarmonic);
 446
 447    double ivn, ivnerror;
 448    Calculate(ivn, ivnerror, fAllStatistics, fPtUComponents, 0, 0);
 449    h->FillIntegratedFlowPOI(ivn, ivnerror);
 450
 451    for(int bin=1; bin<=nBinsPt; ++bin) {
 452       double vn=0;
 453       double evn=0;
 454       if( Calculate(vn, evn, fPtStatistics, fPtUComponents, bin) && evn>0 ) {
 455          h->FillDifferentialFlowPtPOI(bin, vn, evn);
 456       }
 457    }
 458
 459    for(int bin=1; bin<=nBinsEta; ++bin) {
 460       double vn=0;
 461       double evn=0;
 462       if( Calculate(vn, evn, fEtaStatistics, fEtaUComponents, bin) && evn>0 ) {
 463          h->FillDifferentialFlowEtaPOI(bin, vn, evn);
 464       }
 465    }
 466
 467    file->WriteTObject(h,0,"Overwrite");
 468
 469    TH1::AddDirectory(oldAddStatus);    // Restore the automatic storage of histograms to its original status
 470
 471    // Restore AliFlowCommonConstants to make sure that no other analysis are affected
 472    c->SetNbinsPt(oldNbinsPt);
 473    c->SetPtMin(oldPtMin);
 474    c->SetPtMax(oldPtMax);
 475    c->SetNbinsEta(oldNbinsEta);
 476    c->SetEtaMin(oldEtaMin);
 477    c->SetEtaMax(oldEtaMax);
 478 }
 479
 480 TList *AliFlowAnalysisWithMSP::ListHistograms()
 481 {
 482    if( fHistList ) {
 483       fHistList->SetOwner(kFALSE);
 484       delete fHistList;
 485    }
 486    fHistList = new TList();
 487    fHistList->SetOwner(kFALSE);
 488
 489    if(fCommonHist)      fHistList->Add(fCommonHist);        // Standard control histograms, if enabled
 490
 491    //Correlations
 492    if(fQaComponents)    fHistList->Add(fQaComponents);      // Averages of Qa components per event for NUA
 493    if(fQbComponents)    fHistList->Add(fQbComponents);      // Averages of Qb components per event for NUA
 494    if(fQaQb)            fHistList->Add(fQaQb);              // Average of QaQb per event
 495    if(fPtUComponents)   fHistList->Add(fPtUComponents);     // ux and uy per pt bin for NUA
 496    if(fEtaUComponents)  fHistList->Add(fEtaUComponents);    // ux and uy per eta bin for NUA
 497    if(fAllStatistics)   fHistList->Add(fAllStatistics);     // Correlations for uQa uQb and QaQb (integrated)
 498    if(fPtStatistics)    fHistList->Add(fPtStatistics);      // Correlations for uQa uQb and QaQb per pt bin
 499    if(fEtaStatistics)   fHistList->Add(fEtaStatistics);     // Correlations for uQa uQb and QaQb per eta bin
 500
 501    // Result histograms (if calculated)
 502    if(fIntegratedFlow)  fHistList->Add(fIntegratedFlow);    // vn for POI and subevents
 503    if(fDiffFlowPt)      fHistList->Add(fDiffFlowPt);        // vn as function of pt
 504    if(fDiffFlowEta)     fHistList->Add(fDiffFlowEta);       // vn as function of eta
 505    if(fFlags)           fHistList->Add(fFlags);             // Stores fHarmonic and fNUA
 506
 507    return fHistList;
 508 }
 509
 510
 511 void AliFlowAnalysisWithMSP::Print(const Option_t *opt)const
 512 {
 513    if( opt ) std::cout << std::endl;
 514
 515    std::cout << "****************************************************" << std::endl;
 516    std::cout << "    Integrated flow from Modified Scalar Product    " << std::endl;
 517    std::cout << "                                                    " << std::endl;
 518
 519    double vn=0;
 520    double vnerror=0;
 521
 522    std::cout << setprecision(4);
 523    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 0);
 524    std::cout << "v" << fHarmonic << " for POI       : " << setw(11) << vn << " +- " << setw(9) << vnerror << std::endl;
 525    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 1);
 526    std::cout << "v" << fHarmonic << " for subevent A: " << setw(11) << vn << " +- " << setw(9) << vnerror << std::endl;
 527    Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 2);
 528    std::cout << "v" << fHarmonic << " for subevent B: " << setw(11) << vn << " +- " << setw(9) << vnerror << std::endl;
 529    std::cout << std::endl;
 530
 531    std::cout << "NUA terms: " << (fNUA?"(applied)":"(NOT applied)") << std::endl;
 532    std::cout << setprecision(3);
 533    const double ux=fPtUComponents->Average(0);       // Average over all bins
 534    const double eux=TMath::Sqrt(fPtUComponents->Variance(0));
 535    std::cout << "<ux>       " << setw(12) << ux << " +- " << setw(12) << eux << (TMath::Abs(ux)<2*eux?" NOT significant ":" ") << std::endl;
 536    const double ux0eta=fEtaUComponents->Average(fEtaUComponents->FindBin(0.0),0);
 537    const double eux0eta=TMath::Sqrt(fEtaUComponents->Variance(fEtaUComponents->FindBin(0.0),0));
 538    std::cout << "<ux> eta=0 " << setw(12) << ux0eta << " +- " <<  setw(12) << eux0eta << (TMath::Abs(ux0eta)<2*eux0eta?" NOT significant":" ") << std::endl;
 539    const double ux0pt=fPtUComponents->Average(fPtUComponents->FindBin(1.0),0);
 540    const double eux0pt=TMath::Sqrt(fPtUComponents->Variance(fPtUComponents->FindBin(1.0),0));
 541    std::cout << "<ux> pt=1  " << setw(12) << ux0pt << " +- " <<  setw(12) << eux0pt << (TMath::Abs(ux0pt)<2*eux0pt?" NOT significant":" ") << std::endl;
 542
 543    const double uy=fPtUComponents->Average(0);        // Average over all bins
 544    const double euy=TMath::Sqrt(fPtUComponents->Variance(0));
 545    std::cout << "<uy>       " << setw(12) << uy << " +- " << setw(12) << euy << (TMath::Abs(uy)<2*euy?" NOT significant ":" ") << std::endl;;
 546    const double uy0eta=fEtaUComponents->Average(fEtaUComponents->FindBin(0.0),1);
 547    const double euy0eta=TMath::Sqrt(fEtaUComponents->Variance(fEtaUComponents->FindBin(0.0),1));
 548    std::cout << "<uy> eta=0 " << setw(12) << uy0eta << " +- " << setw(12) << euy0eta << (TMath::Abs(uy0eta)<2*euy0eta?" NOT significant ":" ") << std::endl;
 549    const double uy0pt=fPtUComponents->Average(fPtUComponents->FindBin(1.0),1);
 550    const double euy0pt=TMath::Sqrt(fPtUComponents->Variance(fPtUComponents->FindBin(1.0),1));
 551    std::cout << "<uy> pt=1  " << setw(12) << uy0pt << " +- " << setw(12) << euy0pt << (TMath::Abs(uy0pt)<2*euy0pt?" NOT significant ":" ") << std::endl;
 552
 553    const double ax=fQaComponents->Average(0);
 554    const double eax=TMath::Sqrt(fQaComponents->Variance(0));
 555    std::cout << "<QaX>      " << setw(12) << ax << " +- " << setw(12) <<eax << (TMath::Abs(ax)<2*eax?" NOT significant ":" ") << std::endl;
 556    const double ay=fQaComponents->Average(1);
 557    const double eay=TMath::Sqrt(fQaComponents->Variance(1));
 558    std::cout << "<QaY>      " << setw(12) << ay << " +- " << setw(12) << eay << (TMath::Abs(ay)<2*eay?" NOT significant ":" ") << std::endl;
 559    const double bx=fQbComponents->Average(0);
 560    const double ebx=TMath::Sqrt(fQbComponents->Variance(0));
 561    std::cout << "<QbX>      " << setw(12) << bx << " +- " << setw(12) << ebx << (TMath::Abs(bx)<2*ebx?" NOT significant ":" ") << std::endl;
 562    const double by=fQbComponents->Average(1);
 563    const double eby=TMath::Sqrt(fQbComponents->Variance(1));
 564    std::cout << "<QbY>      " << setw(12) << by << " +- " << setw(12) << eby << (TMath::Abs(by)<2*eby?" NOT significant ":" ") << std::endl;
 565    const double ab=fQaQb->Average(0);
 566    const double eab=TMath::Sqrt(fQbComponents->Variance(0));
 567    std::cout << "<QaQb>     " << setw(12) << ab << " +- " << setw(12) << eab << (TMath::Abs(ab)<2*eab?" NOT significant ":" ") << std::endl;
 568    std::cout << std::endl;
 569
 570    std::cout << "Covariance matrix: " << std::endl;
 571    std::cout << "     " << setw(12) << "uQa"                           << setw(12) << "uQb"                           << setw(12) << "QaQb" << std::endl;
 572    std::cout << "uQa  " << setw(12) << fAllStatistics->Covariance(0,0) << std::endl;
 573    std::cout << "uQb  " << setw(12) << fAllStatistics->Covariance(1,0) << setw(12) << fAllStatistics->Covariance(1,1) << std::endl;
 574    std::cout << "QaQb " << setw(12) << fAllStatistics->Covariance(2,0) << setw(12) << fAllStatistics->Covariance(2,1) << setw(12) << fQaQb->Variance(0) << std::endl;
 575    std::cout << "****************************************************" << std::endl;
 576    std::cout << std::endl;
 577 }
 578
 579
 580 AliFlowAnalysisWithMSP &AliFlowAnalysisWithMSP::operator=(const AliFlowAnalysisWithMSP &x)
 581 {
 582    if (&x==this) return *this; //handle self assignmnet
 583    SetNameTitle("MSP","Flow analysis with the Modified Scalar Product method");
 584    delete fQaComponents; fQaComponents=0;
 585    if( x.fQaComponents )   fQaComponents=(AliFlowMSPHistograms *)(x.fQaComponents)->Clone();
 586    delete fQbComponents; fQbComponents=0;
 587    if( x.fQbComponents )   fQbComponents=(AliFlowMSPHistograms *)(x.fQbComponents)->Clone();
 588    delete fQaQb; fQaQb=0;
 589    if( x.fQaQb )           fQaQb=(AliFlowMSPHistograms *)(x.fQaQb)->Clone();
 590    delete fPtUComponents; fPtUComponents=0;
 591    if( fPtUComponents)     fPtUComponents=(AliFlowMSPHistograms *)(x.fPtUComponents)->Clone();
 592    delete fEtaUComponents; fEtaUComponents=0;
 593    if( fEtaUComponents )   fEtaUComponents=(AliFlowMSPHistograms *)(x.fEtaUComponents)->Clone();
 594    delete fAllStatistics; fAllStatistics=0;
 595    if( fAllStatistics )    fAllStatistics=(AliFlowMSPHistograms *)(x.fAllStatistics)->Clone();
 596    delete fPtStatistics; fPtStatistics=0;
 597    if( fPtStatistics )     fPtStatistics=(AliFlowMSPHistograms *)(x.fPtStatistics)->Clone();
 598    delete fEtaStatistics; fEtaStatistics=0;
 599    if( fEtaStatistics )    fEtaStatistics=(AliFlowMSPHistograms *)(x.fEtaStatistics)->Clone();
 600    delete fIntegratedFlow; fIntegratedFlow=0;
 601    if( fIntegratedFlow )   fIntegratedFlow=(TH1D *)(x.fIntegratedFlow)->Clone();
 602    delete fDiffFlowPt; fDiffFlowPt=0;
 603    if( fDiffFlowPt )       fDiffFlowPt=(TH1D *)(x.fDiffFlowPt)->Clone();
 604    delete fDiffFlowEta; fDiffFlowEta=0;
 605    if( fDiffFlowEta )      fDiffFlowEta=(TH1D *)(x.fDiffFlowEta)->Clone();
 606    delete fFlags; fFlags=0;
 607    if( fFlags )            fFlags=(TProfile *)(x.fFlags)->Clone();
 608    delete fCommonHist; fCommonHist=0;
 609    if( fCommonHist ) fCommonHist=new AliFlowCommonHist(*(x.fCommonHist));
 610    return *this;
 611 }
 612
 613 // private functions --------------------------------------------------------------------------------------
 614 bool AliFlowAnalysisWithMSP::Calculate(double &vn, double &vnerror, const AliFlowMSPHistograms *hist, const AliFlowMSPHistograms *components, const int bin, const int poi) const
 615 {
 616    // Get all averages and correlations need for the flow calculation
 617    double uQa=hist->Average(bin,0);                // <<uQa>>
 618    double VuQa=hist->Variance(bin,0);              // Var(<<uQa>>)
 619    double uQb=hist->Average(bin,1);                // <<uQb>>
 620    double VuQb=hist->Variance(bin,1);              // Var(<<uQb>>)
 621    double QaQb=fQaQb->Average(1,0);                // <QaQb> Should not be taken from hist(bin) because there QaQb is entered multiple times: <<QaQb>>!!
 622    double VQaQb=fQaQb->Variance(1,0);              // V(<QaQb>)
 623    double CuQauQb=hist->Covariance(bin,0,1);       // Cov(<<uQa>>,<<uQb>>)
 624    double CuQaQaQb=hist->Covariance(bin,0,2);      // Cov(<<uQa>>,<QaQb>)
 625    double CuQbQaQb=hist->Covariance(bin,1,2);      // Cov(<<uQb>>,<QaQb>)
 626
 627    if( fNUA && components ) {
 628       // Apply NUA correction to QaQb.
 629       double QaX=fQaComponents->Average(0);
 630       double QaY=fQaComponents->Average(1);
 631       double QbX=fQbComponents->Average(0);
 632       double QbY=fQbComponents->Average(1);
 633
 634       QaQb=QaQb-QaX*QbX-QaY*QbY;
 635
 636       // Apply NUA to uQa and uQb (per bin)
 637       double uX=components->Average(bin,0);  // bin 0 is integrated over all bins
 638       double uY=components->Average(bin,1);
 639
 640       uQa = uQa - uX*QaX - uY*QaY;
 641       uQb = uQb - uX*QbX - uY*QbY;
 642       // Error calculation not fully modified: only above terms, the spread in <<u>> and <Qa> and <Qb> are not used
 643       // therefore this should only be applied if significant!
 644       // Check if not fully NUA correcting the error calculation is justified (but this is compatible with the original SP method)!
 645       // In general it is not justified but this can only be checked by splitting the event sample in many subsamples and looking at
 646       // the variance of the result. This may not be feasible if statistics is low
 647    }
 648
 649    // Some sanity checks:
 650    if( uQa*uQb*QaQb <= 0 ) {        // Catch imaginary results
 651       vn=-99;
 652       vnerror=-99;
 653       return false;
 654    }
 655
 656    // Sanity checks passed, calculate, print and store
 657    switch (poi) {
 658    case 1:                                        // Subevent A reference flow
 659       {
 660          if( TMath::Abs(uQb) < 1e-30*TMath::Abs(uQa*QaQb) ) {  // Protect against infinity
 661             vn=0;
 662             vnerror=-1;
 663             return false;
 664          }
 665          double vnA = TMath::Sqrt( uQa*QaQb / uQb ); // vn
 666          double VvnA = QaQb*VuQa/(4*uQa*uQb)         // Variance of vn
 667             + uQa*VQaQb/(4*QaQb*uQb)
 668             + uQa*QaQb*VuQb/(4*TMath::Power(uQb,3))
 669             + CuQaQaQb/(2*uQb)
 670             - QaQb*CuQauQb/(2*TMath::Power(uQb,2))
 671             - uQa*CuQbQaQb/(2*TMath::Power(uQb,2));
 672          vn=vnA;
 673          if( VvnA<0 ) {
 674             vnerror=VvnA;
 675             return false;
 676          }
 677          vnerror=TMath::Sqrt(VvnA);
 678       }
 679       break;
 680    case 2:                                        // Subevent B reference flow
 681       {
 682          if( TMath::Abs(uQa) < 1e-30*TMath::Abs(uQb*QaQb) ) {  // Protect against infinity
 683             vn=0;
 684             vnerror=-1;
 685             return false;
 686          }
 687          double vnB = TMath::Sqrt( uQb*QaQb / uQa );        // vn
 688          double VvnB = uQb*VQaQb/(4*QaQb*uQa)               // Variance of vn
 689             + QaQb*VuQb/(4*uQb*uQa)
 690             + QaQb*uQb*VuQa/(4*TMath::Power(uQa,3))
 691             + CuQbQaQb/(2*uQa)
 692             - uQb*CuQaQaQb/(2*TMath::Power(uQa,2))
 693             - QaQb*CuQauQb/(2*TMath::Power(uQa,2));
 694          vn=vnB;
 695          if( VvnB<0 ) {
 696             vnerror=VvnB;
 697             return false;
 698          }
 699          vnerror=TMath::Sqrt(VvnB);
 700       }
 701       break;
 702    default:                                       // POI flow
 703       {
 704          if( TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQa*uQb) ) {  // Catch infinity
 705             vn=0;
 706             vnerror=-1;
 707             return false;
 708          }
 709          double vnP = TMath::Sqrt( uQa*uQb / QaQb );        // vn
 710          if(   TMath::Abs(uQa*QaQb) < 1e-30*TMath::Abs(uQb*VuQa)
 711             || TMath::Abs(uQb*QaQb) < 1e-30*TMath::Abs(uQa*VuQb)
 712             || TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQa*uQb*VQaQb)
 713             || TMath::Abs(QaQb) < 1e-30*TMath::Abs(CuQauQb)
 714             || TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQb*CuQaQaQb)
 715             || TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQa*CuQbQaQb)
 716             ){
 717                vnerror=-98;
 718                return false;
 719          }
 720          double VvnP = uQb*VuQa/(4*uQa*QaQb)                // Variance of vn
 721             + uQa*VuQb/(4*uQb*QaQb)
 722             + uQa*uQb*VQaQb/(4*TMath::Power(QaQb,3))
 723             + CuQauQb/(2*QaQb)
 724             - uQb*CuQaQaQb/(2*TMath::Power(QaQb,2))
 725             - uQa*CuQbQaQb/(2*TMath::Power(QaQb,2));
 726          vn=TMath::Sign(vnP,uQb);
 727          if( VvnP<0 ) {
 728             vnerror=VvnP;
 729             return false;
 730          }
 731          vnerror=TMath::Sqrt(VvnP);
 732       }
 733    }  // Switch between POI and subevents
 734
 735    return (vnerror>=0);
 736 }
 737
 738 void AliFlowAnalysisWithMSP::ReadHistograms(TDirectory *file)
 739 {
 740    delete fCommonHist;     fCommonHist=0;    // Delete existing histograms
 741    delete fQaComponents;   fQaComponents=0;
 742    delete fQbComponents;   fQbComponents=0;
 743    delete fQaQb;           fQaQb=0;
 744    delete fPtUComponents;  fPtUComponents=0;
 745    delete fEtaUComponents; fEtaUComponents=0;
 746    delete fAllStatistics;  fAllStatistics=0;
 747    delete fPtStatistics;   fPtStatistics=0;
 748    delete fEtaStatistics;  fEtaStatistics=0;
 749    delete fIntegratedFlow; fIntegratedFlow=0;
 750    delete fDiffFlowPt;     fDiffFlowPt=0;
 751    delete fDiffFlowEta;    fDiffFlowEta=0;
 752    delete fFlags;          fFlags=0;
 753    if( fHistList ) {
 754       fHistList->SetOwner(kFALSE);
 755       delete fHistList;
 756       fHistList=0;
 757    }
 758
 759    file->GetObject("QaComponents",fQaComponents);
 760    file->GetObject("QbComponents",fQbComponents);
 761    file->GetObject("QaQb",fQaQb);
 762    file->GetObject("PtUComponents",fPtUComponents);
 763    file->GetObject("EtaUComponents",fEtaUComponents);
 764    file->GetObject("AllStatistics",fAllStatistics);
 765    file->GetObject("PtStatistics",fPtStatistics);
 766    file->GetObject("EtaStatistics",fEtaStatistics);
 767    if( !fQaComponents || !fQbComponents || !fQaQb || !fPtUComponents || !fEtaUComponents || !fAllStatistics || !fPtStatistics || !fEtaStatistics ) {
 768       std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms() : One or more histograms were not read correctly from " << file->GetPath() << std::endl;
 769    }
 770
 771    file->GetObject("Flags",fFlags);          // Flags are required
 772
 773    if( !fFlags ){
 774       std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *) : Flags histogram not found, using defaults" << std::endl;
 775       fHarmonic=2;
 776       fNUA=false;
 777    }else{
 778       fHarmonic=(UInt_t)(fFlags->GetBinContent(1));
 779       double harmonicSpread=fFlags->GetBinError(1);
 780       if( harmonicSpread!=0 ) {
 781          std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *) :These histograms seem to be merged from analysis with different Harmonics. Results removed!" << std::endl;
 782          delete fIntegratedFlow; fIntegratedFlow=0;
 783          delete fDiffFlowPt;     fDiffFlowPt=0;
 784          delete fDiffFlowEta;    fDiffFlowEta=0;
 785       }
 786       fNUA=fFlags->GetBinContent(2);   // Mixing NUA does not matter since it needs to be recalculated anyway
 787       double nuaSpread=fFlags->GetBinError(2);
 788       if( nuaSpread!=0 ) {
 789          std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *) :These histograms seem to be merged from analysis with different NUA corrections. Results removed" << std::endl;
 790          delete fIntegratedFlow; fIntegratedFlow=0;
 791          delete fDiffFlowPt;     fDiffFlowPt=0;
 792          delete fDiffFlowEta;    fDiffFlowEta=0;
 793       }
 794    }
 795
 796    // Optional histograms, may return a zero pointer
 797    file->GetObject("AliFlowCommonHist_MSP",fCommonHist);             // The AliFlowCommonHist is optional
 798    file->GetObject("IntegratedFlow",fIntegratedFlow);                // Results are optional
 799    file->GetObject("DiffFlowPt",fDiffFlowPt);
 800    file->GetObject("DiffFlowEta",fDiffFlowEta);
 801
 802    fBookCommonHistograms=(fCommonHist!=0);
 803 }
 804
 805
 806 void AliFlowAnalysisWithMSP::ReadHistograms(TList *list)
 807 {
 808    if( !list ) return;
 809    delete fCommonHist;     fCommonHist=0;    // Delete existing histograms if any
 810    delete fQaComponents;   fQaComponents=0;
 811    delete fQbComponents;   fQbComponents=0;
 812    delete fQaQb;           fQaQb=0;
 813    delete fPtUComponents;  fPtUComponents=0;
 814    delete fEtaUComponents; fEtaUComponents=0;
 815    delete fAllStatistics;  fAllStatistics=0;
 816    delete fPtStatistics;   fPtStatistics=0;
 817    delete fEtaStatistics;  fEtaStatistics=0;
 818    delete fIntegratedFlow; fIntegratedFlow=0;
 819    delete fDiffFlowPt;     fDiffFlowPt=0;
 820    delete fDiffFlowEta;    fDiffFlowEta=0;
 821    delete fFlags;          fFlags=0;
 822    if( fHistList ) {
 823       fHistList->SetOwner(kFALSE);
 824       delete fHistList;
 825       fHistList=0;
 826    }
 827
 828    fQaComponents = static_cast<AliFlowMSPHistograms *>(list->FindObject("QaComponents"));
 829    fQbComponents = static_cast<AliFlowMSPHistograms *>(list->FindObject("QbComponents"));
 830    fQaQb = static_cast<AliFlowMSPHistograms *>(list->FindObject("QaQb"));
 831    fPtUComponents = static_cast<AliFlowMSPHistograms *>(list->FindObject("PtUComponents"));
 832    fEtaUComponents = static_cast<AliFlowMSPHistograms *>(list->FindObject("EtaUComponents"));
 833    fAllStatistics = static_cast<AliFlowMSPHistograms *>(list->FindObject("AllStatistics"));
 834    fPtStatistics = static_cast<AliFlowMSPHistograms *>(list->FindObject("PtStatistics"));
 835    fEtaStatistics = static_cast<AliFlowMSPHistograms *>(list->FindObject("EtaStatistics"));
 836    if( !fQaComponents || !fQbComponents || !fQaQb || !fPtUComponents || !fEtaUComponents || !fAllStatistics || !fPtStatistics || !fEtaStatistics ) {
 837       std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms(Tlist *) : One or more histograms were not read correctly from TList" << std::endl;
 838    }
 839
 840    fFlags = static_cast<TProfile *>(list->FindObject("Flags"));          // Flags are required
 841
 842    if( !fFlags ){
 843       std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms(TList *) : Flags histogram not found, using defaults" << std::endl;
 844       fHarmonic=2;
 845       fNUA=false;
 846    }else{
 847       fHarmonic=(UInt_t)(fFlags->GetBinContent(1));
 848       double harmonicSpread=fFlags->GetBinError(1);
 849       if( harmonicSpread!=0 ) {
 850          std::cerr << "These histograms seem to be merged from analysis with different Harmonics. Results removed!" << std::endl;
 851          delete fIntegratedFlow; fIntegratedFlow=0;
 852          delete fDiffFlowPt;     fDiffFlowPt=0;
 853          delete fDiffFlowEta;    fDiffFlowEta=0;
 854       }
 855       fNUA=fFlags->GetBinContent(2);   // Mixing NUA does not matter since it needs to be recalculated anyway
 856       double nuaSpread=fFlags->GetBinError(2);
 857       if( nuaSpread!=0 ) {
 858          std::cerr << "These histograms seem to be merged from analysis with different NUA corrections. Results removed" << std::endl;
 859          delete fIntegratedFlow; fIntegratedFlow=0;
 860          delete fDiffFlowPt;     fDiffFlowPt=0;
 861          delete fDiffFlowEta;    fDiffFlowEta=0;
 862       }
 863    }
 864
 865    // Optional histograms, may return a zero pointer
 866    fCommonHist = static_cast<AliFlowCommonHist *>(list->FindObject("AliFlowCommonHist_MSP"));             // The AliFlowCommonHist is optional
 867    fIntegratedFlow = static_cast<TH1D *>(list->FindObject("IntegratedFlow"));                // Results are optional
 868    fDiffFlowPt = static_cast<TH1D *>(list->FindObject("DiffFlowPt"));
 869    fDiffFlowEta = static_cast<TH1D *>(list->FindObject("DiffFlowEta"));
 870
 871    fBookCommonHistograms=(fCommonHist!=0);
 872 }