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