PWG/FLOW/Base/AliFlowAnalysisWithFittingQDistribution.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 /********************************
  17  * estimating reference flow by *
  18  *   fitting q-distribution     *
  19  *                              *
  20  * author: Ante Bilandzic       *
  21  *       (abilandzic@gmail.com) *
  22  *                              *
  23  *  based on the macro written  *
  24  *     by Sergei Voloshin       *
  25  *******************************/
  26
  27 #define AliFlowAnalysisWithFittingQDistribution_cxx
  28
  29 #include "Riostream.h"
  30 #include "AliFlowCommonConstants.h"
  31 #include "AliFlowCommonHist.h"
  32 #include "AliFlowCommonHistResults.h"
  33 #include "TChain.h"
  34 #include "TFile.h"
  35 #include "TList.h"
  36 #include "TF1.h"
  37 #include "TH2D.h"
  38 #include "TParticle.h"
  39 #include "TProfile.h"
  40 #include "AliFlowEventSimple.h"
  41 #include "AliFlowTrackSimple.h"
  42 #include "AliFlowAnalysisWithFittingQDistribution.h"
  43
  44 class TH1;
  45 class TGraph;
  46 class TPave;
  47 class TLatex;
  48 class TMarker;
  49 class TObjArray;
  50 class TList;
  51 class TSystem;
  52 class TROOT;
  53 class AliFlowVector;
  54 class TVector;
  55
  56 //================================================================================================================
  57
  58 using std::endl;
  59 using std::cout;
  60 ClassImp(AliFlowAnalysisWithFittingQDistribution)
  61
  62 AliFlowAnalysisWithFittingQDistribution::AliFlowAnalysisWithFittingQDistribution():
  63  fHistList(NULL),
  64  fBookOnlyBasicCCH(kTRUE),
  65  fCommonHists(NULL),
  66  fCommonHistsResults(NULL),
  67  fnBinsPhi(0),
  68  fPhiMin(0),
  69  fPhiMax(0),
  70  fPhiBinWidth(0),
  71  fnBinsPt(0),
  72  fPtMin(0),
  73  fPtMax(0),
  74  fPtBinWidth(0),
  75  fnBinsEta(0),
  76  fEtaMin(0),
  77  fEtaMax(0),
  78  fEtaBinWidth(0),
  79  fHarmonic(2),
  80  fAnalysisLabel(NULL),
  81  fMultiplicityIs(AliFlowCommonConstants::kRP),
  82  fWeightsList(NULL),
  83  fUsePhiWeights(kFALSE),
  84  fUsePtWeights(kFALSE),
  85  fUseEtaWeights(kFALSE),
  86  fUseParticleWeights(NULL),
  87  fPhiWeights(NULL),
  88  fPtWeights(NULL),
  89  fEtaWeights(NULL),
  90  fSumOfParticleWeights(NULL),
  91  fqDistribution(NULL),
  92  fqMin(0.),
  93  fqMax(100.),
  94  fqNbins(10000),
  95  fStoreqDistributionVsMult(kFALSE),
  96  fqDistributionVsMult(NULL),
  97  fMinMult(0.),
  98  fMaxMult(10000.),
  99  fnBinsMult(1000),
 100  fFittingParameters(NULL),
 101  fTreshold(5.),
 102  fvStart(0.05),
 103  fvMin(0.0),
 104  fvMax(0.25),
 105  fSigma2Start(0.75),
 106  fSigma2Min(0.5),
 107  fSigma2Max(2.5),
 108  fFinalResultIsFromSigma2Fitted(kTRUE),
 109  fPrintOnTheScreen(kTRUE),
 110  fDoFit(kTRUE)
 111  {
 112   // constructor
 113
 114   // base list to hold all output objects:
 115   fHistList = new TList();
 116   fHistList->SetName("cobjFQD");
 117   fHistList->SetOwner(kTRUE);
 118
 119   // analysis label;
 120   fAnalysisLabel = new TString();
 121
 122   // list to hold histograms with phi, pt and eta weights:
 123   fWeightsList = new TList();
 124
 125   // initialize all arrays:
 126   this->InitializeArrays();
 127
 128  } // end of constructor
 129
 130 //================================================================================================================
 131
 132 AliFlowAnalysisWithFittingQDistribution::~AliFlowAnalysisWithFittingQDistribution()
 133 {
 134  // destructor
 135
 136  delete fHistList;
 137
 138 } // end of destructor
 139
 140 //================================================================================================================
 141
 142
 143 void AliFlowAnalysisWithFittingQDistribution::Init()
 144 {
 145  // Access constants and book everything.
 146
 147  //save old value and prevent histograms from being added to directory
 148  //to avoid name clashes in case multiple analaysis objects are used
 149  //in an analysis
 150  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 151  TH1::AddDirectory(kFALSE);
 152
 153  // access constants:
 154  this->AccessConstants();
 155
 156  // booking:
 157  this->BookCommonHistograms();
 158  this->BookAndFillWeightsHistograms();
 159  this->BookEverythingForDistributions();
 160
 161  // store fitting parameters:
 162  this->StoreFittingParameters();
 163
 164  // nest lists:
 165  fWeightsList->SetName("Weights");
 166  fWeightsList->SetOwner(kTRUE);
 167  fHistList->Add(fWeightsList);
 168
 169  // set harmonic in common control histograms (to be improved (should I do this somewhere else?)):
 170  (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic);
 171
 172  //restore old status
 173  TH1::AddDirectory(oldHistAddStatus);
 174 } // end of void AliFlowAnalysisWithFittingQDistribution::Init()
 175
 176 //================================================================================================================
 177
 178 void AliFlowAnalysisWithFittingQDistribution::Make(AliFlowEventSimple* anEvent)
 179 {
 180  // Loop over data only in this method.
 181
 182  // a) Check all pointers used in this method;
 183  // b) Fill the common control histograms;
 184  // c) Loop over data and calculate Q-vector and sum of particle weights;
 185  // d) Fill the histogram for q-distribution and sum of particle weights.
 186
 187  // a) Check all pointers used in this method:
 188  this->CheckPointersUsedInMake();
 189
 190  // b) fill the common control histograms:
 191  fCommonHists->FillControlHistograms(anEvent);
 192
 193  // c) Loop over data and fill histogram for q-distribution:
 194  Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
 195  Double_t dPt  = 0.; // transverse momentum
 196  Double_t dEta = 0.; // pseudorapidity
 197  Double_t wPhi = 1.; // phi weight
 198  Double_t wPt  = 1.; // pt weight
 199  Double_t wEta = 1.; // eta weight
 200  Double_t dReQ = 0.; // real part of Q-vector
 201  Double_t dImQ = 0.; // imaginary part of Q-vector
 202  Int_t n = fHarmonic; // shortcut for the harmonic
 203  Double_t dSumOfParticleWeights = 0.; // when particle weights are not used dSumOfParticleWeights is equal to multiplicity
 204  AliFlowTrackSimple *aftsTrack = NULL;
 205  Int_t nPrim = anEvent->NumberOfTracks(); // nPrim = total number of primary tracks
 206
 207  // Start loop over particles:
 208  for(Int_t i=0;i<nPrim;i++)
 209  {
 210   aftsTrack=anEvent->GetTrack(i);
 211   if(aftsTrack)
 212   {
 213    if(!(aftsTrack->InRPSelection())){continue;} // consider only tracks which are RPs
 214    dPhi = aftsTrack->Phi();
 215    dPt  = aftsTrack->Pt();
 216    dEta = aftsTrack->Eta();
 217    if(fUsePhiWeights && fPhiWeights && fnBinsPhi) // determine phi weight for this particle:
 218    {
 219     wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
 220    }
 221    if(fUsePtWeights && fPtWeights && fPtBinWidth) // determine pt weight for this particle:
 222    {
 223     wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth)));
 224    }
 225    if(fUseEtaWeights && fEtaWeights && fEtaBinWidth) // determine eta weight for this particle:
 226    {
 227     wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth)));
 228    }
 229    // Calculate real and imaginary part of Q-vector and sum of particle weights for this event:
 230    // Q-vector:
 231    dReQ += wPhi*wPt*wEta*TMath::Cos(n*dPhi);
 232    dImQ += wPhi*wPt*wEta*TMath::Sin(n*dPhi);
 233    // sum of particle weights:
 234    dSumOfParticleWeights += wPhi*wPt*wEta; // when particle weights are not used this sum gives # of RPs, i.e. multiplicity
 235   } // end of if(aftsTrack)
 236  } // end of for(Int_t i=0;i<nPrim;i++)
 237
 238  // d) Fill the histogram for q-distribution and sum of particle weights:
 239  Double_t q = 0.; // q = Q\sqrt{sum of particle weights}
 240  if(dSumOfParticleWeights > 1.)
 241  {
 242   q = pow(dReQ*dReQ+dImQ*dImQ,0.5)/pow(dSumOfParticleWeights,0.5);
 243   fqDistribution->Fill(q,1.);
 244   fSumOfParticleWeights->Fill(dSumOfParticleWeights,1.);
 245   if(fStoreqDistributionVsMult)
 246   {
 247    // Multiplicity bin of an event (relevant for all histos vs M):
 248    Double_t dMultiplicityBin = 0.;
 249    if(fMultiplicityIs==AliFlowCommonConstants::kRP)
 250    {
 251     Int_t nNumberOfRPsEBE = anEvent->GetNumberOfRPs(); // number of RPs (i.e. number of reference particles)
 252     dMultiplicityBin = nNumberOfRPsEBE+0.5;
 253    } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
 254      {
 255       Int_t nReferenceMultiplicityEBE = anEvent->GetReferenceMultiplicity(); // reference multiplicity for current event
 256       dMultiplicityBin = nReferenceMultiplicityEBE+0.5;
 257      } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
 258        {
 259         Int_t nNumberOfPOIsEBE = anEvent->GetNumberOfPOIs(); // number of POIs (i.e. number of particles of interest)
 260         dMultiplicityBin = nNumberOfPOIsEBE+0.5;
 261        }
 262    // Fill qDist vs M:
 263    fqDistributionVsMult->Fill(q,dMultiplicityBin);
 264   } // end of if(fStoreqDistributionVsMult)
 265  } // end of if(dSumOfParticleWeights > 1.)
 266
 267 } // end of Make()
 268
 269 //================================================================================================================
 270
 271 void AliFlowAnalysisWithFittingQDistribution::Finish(Bool_t doFit)
 272 {
 273  // Calculate the final results.
 274
 275  // a) Check all pointers used in this method;
 276  // b) Acces common constants;
 277  // c) Access fitting paremeters;
 278  // d) Do the final fit of q-distribution;
 279  // e) Fill common hist results;
 280  // f) Print on the screen the final results.
 281
 282  // a) Check all pointers used in this method:
 283  this->CheckPointersUsedInFinish();
 284
 285  // b) Access common constants:
 286  this->AccessConstants();
 287
 288  // c) Access fitting paremeters:
 289  this->AccessFittingParameters();
 290
 291  // d) Do the final fit of q-distribution:
 292  if(doFit && fDoFit)
 293  {
 294   this->DoFit(kFALSE); // sigma^2 not fitted (fixed to 0.5)
 295   this->DoFit(kTRUE); // sigma^2 fitted
 296  }
 297
 298  // e) Fill common hist results (by default fill with results obtained with sigma^2 fitted,
 299  //    alternatively use a setter SetFinalResultIsFromSigma2Fitted(kFALSE)):
 300  this->FillCommonHistResults(fFinalResultIsFromSigma2Fitted);
 301
 302  // f) Print on the screen the final results:
 303  if(fPrintOnTheScreen) this->PrintOnTheScreen();
 304
 305 } // end of void AliFlowAnalysisWithFittingQDistribution::Finish(Bool_t doFit)
 306
 307 //================================================================================================================
 308
 309 void AliFlowAnalysisWithFittingQDistribution::CheckPointersUsedInMake()
 310 {
 311  // Check all pointers used in method Make().
 312
 313  if(!fCommonHists)
 314  {
 315   cout<<endl;
 316   cout<<" WARNING (FQD::Make()): fCommonHists is NULL !!!!"<<endl;
 317   cout<<endl;
 318   exit(0);
 319  }
 320  if(!fqDistribution)
 321  {
 322   cout<<endl;
 323   cout<<" WARNING (FQD::Make()): fqDistribution is NULL !!!!"<<endl;
 324   cout<<endl;
 325   exit(0);
 326  }
 327  if(!fSumOfParticleWeights)
 328  {
 329   cout<<endl;
 330   cout<<" WARNING (FQD::Make()): fSumOfParticleWeights is NULL !!!!"<<endl;
 331   cout<<endl;
 332   exit(0);
 333  }
 334  if(fUsePhiWeights && !fPhiWeights)
 335  {
 336   cout<<endl;
 337   cout<<" WARNING (FQD::Make()): fPhiWeights is NULL !!!!"<<endl;
 338   cout<<endl;
 339   exit(0);
 340  }
 341  if(fUsePtWeights && !fPtWeights)
 342  {
 343   cout<<endl;
 344   cout<<" WARNING (FQD::Make()): fPtWeights is NULL !!!!"<<endl;
 345   cout<<endl;
 346   exit(0);
 347  }
 348  if(fUseEtaWeights && !fEtaWeights)
 349  {
 350   cout<<endl;
 351   cout<<" WARNING (FQD::Make()): fEtaWeights is NULL !!!!"<<endl;
 352   cout<<endl;
 353   exit(0);
 354  }
 355
 356 } // end of void AliFlowAnalysisWithFittingQDistribution::CheckPointersUsedInMake()
 357
 358 //================================================================================================================
 359
 360 void AliFlowAnalysisWithFittingQDistribution::CheckPointersUsedInFinish()
 361 {
 362  // Check all pointers used in method Finish().
 363
 364  if(!fFittingParameters)
 365  {
 366   cout<<endl;
 367   cout<<" WARNING (FQD::Finish()): fFittingParameters is NULL !!!!"<<endl;
 368   cout<<endl;
 369   exit(0);
 370  }
 371  if(!fqDistribution)
 372  {
 373   cout<<endl;
 374   cout<<" WARNING (FQD::Finish()): fqDistribution is NULL !!!!"<<endl;
 375   cout<<endl;
 376   exit(0);
 377  }
 378  if(!fSumOfParticleWeights)
 379  {
 380   cout<<endl;
 381   cout<<" WARNING (FQD::Finish()): fSumOfParticleWeights is NULL !!!!"<<endl;
 382   cout<<endl;
 383   exit(0);
 384  }
 385  for(Int_t s2F=0;s2F<2;s2F++)
 386  {
 387   if(!fIntFlow[s2F])
 388   {
 389    cout<<endl;
 390    cout<<" WARNING (FQD::Finish()): "<<Form("fIntFlow[%d] is NULL !!!!",s2F)<<endl;
 391    cout<<endl;
 392    exit(0);
 393   }
 394   if(!fSigma2[s2F])
 395   {
 396    cout<<endl;
 397    cout<<" WARNING (FQD::Finish()): "<<Form("fSigma2[%d] is NULL !!!!",s2F)<<endl;
 398    cout<<endl;
 399    exit(0);
 400   }
 401   if(!fChi2[s2F])
 402   {
 403    cout<<endl;
 404    cout<<" WARNING (FQD::Finish()): "<<Form("fChi2[%d] is NULL !!!!",s2F)<<endl;
 405    cout<<endl;
 406    exit(0);
 407   }
 408  } // end of for(Int_t s2F=0;s2F<2;s2F++)
 409  if(!fCommonHistsResults)
 410  {
 411   cout<<endl;
 412   cout<<"WARNING (FQD::Finish()): fCommonHistsResults is NULL !!!!"<<endl;
 413   cout<<endl;
 414   exit(0);
 415  }
 416  if(!fCommonHists)
 417  {
 418   cout<<endl;
 419   cout<<"WARNING (FQD::Finish()): fCommonHists is NULL !!!!"<<endl;
 420   cout<<endl;
 421   exit(0);
 422  }
 423
 424 } // end of void AliFlowAnalysisWithFittingQDistribution::CheckPointersUsedInFinish()
 425
 426 //================================================================================================================
 427
 428 void AliFlowAnalysisWithFittingQDistribution::GetOutputHistograms(TList *outputListHistos)
 429 {
 430  // Get pointers to all output histograms (called before Finish()).
 431
 432  if(outputListHistos)
 433  {
 434   // 1.) common control histograms and common histograms for final results:
 435   TString commonHistName = "AliFlowCommonHistFQD";
 436   commonHistName += fAnalysisLabel->Data();
 437   AliFlowCommonHist *commonHist = dynamic_cast<AliFlowCommonHist*>(outputListHistos->FindObject(commonHistName.Data()));
 438   if(commonHist) this->SetCommonHists(commonHist);
 439
 440   TString commonHistResName = "AliFlowCommonHistResultsFQD";
 441   commonHistResName += fAnalysisLabel->Data();
 442   AliFlowCommonHistResults *commonHistRes = dynamic_cast<AliFlowCommonHistResults*>
 443                                             (outputListHistos->FindObject(commonHistResName.Data()));
 444   if(commonHistRes) this->SetCommonHistsResults(commonHistRes);
 445
 446   // 2.) weights:
 447   TList *weightsList = dynamic_cast<TList*>(outputListHistos->FindObject("Weights"));
 448   if(weightsList){this->SetWeightsList(weightsList);}
 449   Bool_t bUsePhiWeights = kFALSE;
 450   Bool_t bUsePtWeights = kFALSE;
 451   Bool_t bUseEtaWeights = kFALSE;
 452   TString fUseParticleWeightsName = "fUseParticleWeightsFQD";
 453   fUseParticleWeightsName += fAnalysisLabel->Data();
 454   TProfile *useParticleWeights = NULL;
 455   if(weightsList)
 456   {
 457    useParticleWeights = dynamic_cast<TProfile*>(weightsList->FindObject(fUseParticleWeightsName.Data()));
 458   }
 459   if(useParticleWeights)
 460   {
 461    this->SetUseParticleWeights(useParticleWeights);
 462    bUsePhiWeights = (Int_t)useParticleWeights->GetBinContent(1);
 463    this->SetUsePhiWeights(bUsePhiWeights);
 464    bUsePtWeights = (Int_t)useParticleWeights->GetBinContent(2);
 465    this->SetUsePtWeights(bUsePtWeights);
 466    bUseEtaWeights = (Int_t)useParticleWeights->GetBinContent(3);
 467    this->SetUseEtaWeights(bUseEtaWeights);
 468   }
 469
 470   // 3.) distributions and 4.) final results of fitting:
 471   TString sigmaFlag[2] = {"#sigma^{2} not fitted","#sigma^{2} fitted"};
 472   // q-distribution:
 473   TString qDistributionName = "fqDistribution";
 474   qDistributionName += fAnalysisLabel->Data();
 475   // sum of particle weights:
 476   TString sumOfParticleWeightsName = "fSumOfParticleWeights";
 477   sumOfParticleWeightsName += fAnalysisLabel->Data();
 478   // final results for reference flow:
 479   TString intFlowName = "fIntFlowFQD";
 480   intFlowName += fAnalysisLabel->Data();
 481   // sigma^2:
 482   TString sigma2Name = "fSigma2";
 483   sigma2Name += fAnalysisLabel->Data();
 484   // chi^2:
 485   TString chi2Name = "fChi2";
 486   chi2Name += fAnalysisLabel->Data();
 487   // fitting function:
 488   TString fittingFunctionName = "fFittingFunction";
 489   fittingFunctionName += fAnalysisLabel->Data();
 490
 491   TH1D *qDistribution = NULL;
 492   TH1D *sumOfParticleWeights = NULL;
 493   TH1D *intFlow[2] = {NULL};
 494   TH1D *sigma2[2] = {NULL};
 495   TH1D *chi2[2] = {NULL};
 496   TF1 *fittingFunction[2] = {NULL};
 497
 498   // q-distribution:
 499   qDistribution = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s",qDistributionName.Data())));
 500   if(qDistribution)
 501   {
 502    this->SetqDistribution(qDistribution);
 503   } else
 504     {
 505      cout<<"WARNING: qDistribution is NULL in AFAWFQD::GOH() !!!!"<<endl;
 506     }
 507   // sum of particle weights:
 508   sumOfParticleWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s",sumOfParticleWeightsName.Data())));
 509   if(sumOfParticleWeights)
 510   {
 511    this->SetSumOfParticleWeights(sumOfParticleWeights);
 512   } else
 513     {
 514      cout<<"WARNING: sumOfParticleWeights is NULL in AFAWFQD::GOH() !!!!"<<endl;
 515     }
 516   // final results:
 517   for(Int_t f=0;f<2;f++)
 518   {
 519    // final results for reference flow:
 520    intFlow[f] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s, %s",intFlowName.Data(),sigmaFlag[f].Data())));
 521    if(intFlow[f])
 522    {
 523     this->SetIntFlow(intFlow[f],f);
 524    } else
 525      {
 526       cout<<"WARNING: intFlow[f] is NULL in AFAWFQD::GOH() !!!!"<<endl;
 527       cout<<"f  = "<<f<<endl;
 528      }
 529    // sigma^2:
 530    sigma2[f] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s, %s",sigma2Name.Data(),sigmaFlag[f].Data())));
 531    if(sigma2[f])
 532    {
 533     this->SetSigma2(sigma2[f],f);
 534    } else
 535      {
 536       cout<<"WARNING: sigma2[f] is NULL in AFAWFQD::GOH() !!!!"<<endl;
 537       cout<<"f  = "<<f<<endl;
 538      }
 539    // chi^2:
 540    chi2[f] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s, %s",chi2Name.Data(),sigmaFlag[f].Data())));
 541    if(chi2[f])
 542    {
 543     this->SetChi2(chi2[f],f);
 544    } else
 545      {
 546       cout<<"WARNING: chi2[f] is NULL in AFAWFQD::GOH() !!!!"<<endl;
 547       cout<<"f  = "<<f<<endl;
 548      }
 549    // fitting functions:
 550    fittingFunction[f] = dynamic_cast<TF1*>(outputListHistos->FindObject(Form("%s, %s",fittingFunctionName.Data(),sigmaFlag[f].Data())));
 551    if(fittingFunction[f])
 552    {
 553     this->SetFittingFunction(fittingFunction[f],f);
 554    } else
 555      {
 556       cout<<"WARNING: fittingFunction[f] is NULL in AFAWFQD::GOH() !!!!"<<endl;
 557       cout<<"f  = "<<f<<endl;
 558      }
 559   } // end of for(Int_t f=0;f<2;f++)
 560
 561   // 5.) fitting parameters:
 562   // q-distribution:
 563   TString fittingParametersName = "fFittingParameters";
 564   fittingParametersName += fAnalysisLabel->Data();
 565   TProfile *fittingParameters = NULL;
 566   fittingParameters = dynamic_cast<TProfile*>(outputListHistos->FindObject(fittingParametersName.Data()));
 567   if(fittingParameters)
 568   {
 569    this->SetFittingParameters(fittingParameters);
 570   } else
 571     {
 572      cout<<"WARNING:fittingParameters is NULL in AFAWFQD::GOH() !!!!"<<endl;
 573     }
 574
 575  } else // to if(outputListHistos)
 576    {
 577     cout<<"WARNING: outputListHistos is NULL in AFAWFQD::GOH() !!!!"<<endl;
 578     exit(0);
 579    }
 580
 581 } // end of void AliFlowAnalysisWithFittingQDistribution::GetOutputHistograms(TList *outputListHistos)
 582
 583 //================================================================================================================
 584
 585 void AliFlowAnalysisWithFittingQDistribution::WriteHistograms(TString* outputFileName)
 586 {
 587  //store the final results in output .root file
 588  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
 589  //output->WriteObject(fHistList, "cobjFQD","SingleKey");
 590  fHistList->SetName("cobjFQD");
 591  fHistList->SetOwner(kTRUE);
 592  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 593  delete output;
 594 }
 595
 596 //================================================================================================================
 597
 598 void AliFlowAnalysisWithFittingQDistribution::WriteHistograms(TString outputFileName)
 599 {
 600  //store the final results in output .root file
 601  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 602  //output->WriteObject(fHistList, "cobjFQD","SingleKey");
 603  fHistList->SetName("cobjFQD");
 604  fHistList->SetOwner(kTRUE);
 605  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 606  delete output;
 607 }
 608
 609 //================================================================================================================
 610
 611 void AliFlowAnalysisWithFittingQDistribution::WriteHistograms(TDirectoryFile *outputFileName)
 612 {
 613  //store the final results in output .root file
 614  fHistList->SetName("cobjFQD");
 615  fHistList->SetOwner(kTRUE);
 616  outputFileName->Add(fHistList);
 617  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
 618 }
 619
 620 //================================================================================================================
 621
 622 void AliFlowAnalysisWithFittingQDistribution::InitializeArrays()
 623 {
 624  // Initialize all arrays.
 625
 626  for(Int_t s2F=0;s2F<2;s2F++) // sigma^2 not fitted (0) or fitted (1)
 627  {
 628   fIntFlow[s2F] = NULL;
 629   fSigma2[s2F] = NULL;
 630   fChi2[s2F] = NULL;
 631   fFittingFunction[s2F] = NULL;
 632  }
 633
 634 } // end of void AliFlowAnalysisWithFittingQDistribution::InitializeArrays()
 635
 636 //================================================================================================================
 637
 638 void AliFlowAnalysisWithFittingQDistribution::BookCommonHistograms()
 639 {
 640  // Book common histograms.
 641
 642  // common control histogram:
 643  TString commonHistName = "AliFlowCommonHistFQD";
 644  commonHistName += fAnalysisLabel->Data();
 645  fCommonHists = new AliFlowCommonHist(commonHistName.Data(),commonHistName.Data(),fBookOnlyBasicCCH);
 646  fHistList->Add(fCommonHists);
 647
 648  // common histograms for final results:
 649  TString commonHistResName = "AliFlowCommonHistResultsFQD";
 650  commonHistResName += fAnalysisLabel->Data();
 651  fCommonHistsResults = new AliFlowCommonHistResults(commonHistResName.Data(),"",fHarmonic);
 652  fHistList->Add(fCommonHistsResults);
 653
 654 } // end of void AliFlowAnalysisWithFittingQDistribution::BookCommonHistograms()
 655
 656 //================================================================================================================
 657
 658 void AliFlowAnalysisWithFittingQDistribution::BookAndFillWeightsHistograms()
 659 {
 660  // Book and fill histograms which hold phi, pt and eta weights.
 661
 662  if(!fWeightsList)
 663  {
 664   cout<<"WARNING: fWeightsList is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
 665   exit(0);
 666  }
 667
 668  TString fUseParticleWeightsName = "fUseParticleWeightsFQD";
 669  fUseParticleWeightsName += fAnalysisLabel->Data();
 670  fUseParticleWeights = new TProfile(fUseParticleWeightsName.Data(),"0 = particle weight not used, 1 = particle weight used ",3,0,3);
 671  fUseParticleWeights->SetLabelSize(0.08);
 672  (fUseParticleWeights->GetXaxis())->SetBinLabel(1,"w_{#phi}");
 673  (fUseParticleWeights->GetXaxis())->SetBinLabel(2,"w_{p_{T}}");
 674  (fUseParticleWeights->GetXaxis())->SetBinLabel(3,"w_{#eta}");
 675  fUseParticleWeights->Fill(0.5,(Int_t)fUsePhiWeights);
 676  fUseParticleWeights->Fill(1.5,(Int_t)fUsePtWeights);
 677  fUseParticleWeights->Fill(2.5,(Int_t)fUseEtaWeights);
 678  fWeightsList->Add(fUseParticleWeights);
 679  // phi weights:
 680  if(fUsePhiWeights)
 681  {
 682   if(fWeightsList->FindObject("phi_weights"))
 683   {
 684    fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
 685    if(!fPhiWeights){printf("\n WARNING (FQD): !fPhiWeights !!!!\n");exit(0);}
 686    if(TMath::Abs(fPhiWeights->GetBinWidth(1)-fPhiBinWidth)>pow(10.,-6.))
 687    {
 688     cout<<endl;
 689     cout<<"WARNING (FQD): Inconsistent binning in histograms for phi-weights throughout the code."<<endl;
 690     cout<<endl;
 691     //exit(0);
 692    }
 693   } else
 694     {
 695      cout<<"WARNING: fWeightsList->FindObject(\"phi_weights\") is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
 696      exit(0);
 697     }
 698  } // end of if(fUsePhiWeights)
 699  // pt weights:
 700  if(fUsePtWeights)
 701  {
 702   if(fWeightsList->FindObject("pt_weights"))
 703   {
 704    fPtWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("pt_weights"));
 705    if(!fPtWeights){printf("\n WARNING (FQD): !fPtWeights !!!!\n");exit(0);}
 706    if(TMath::Abs(fPtWeights->GetBinWidth(1)-fPtBinWidth)>pow(10.,-6.))
 707    {
 708     cout<<endl;
 709     cout<<"WARNING (FQD): Inconsistent binning in histograms for pt-weights throughout the code."<<endl;
 710     cout<<endl;
 711     //exit(0);
 712    }
 713   } else
 714     {
 715      cout<<"WARNING: fWeightsList->FindObject(\"pt_weights\") is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
 716      exit(0);
 717     }
 718  } // end of if(fUsePtWeights)
 719  // eta weights:
 720  if(fUseEtaWeights)
 721  {
 722   if(fWeightsList->FindObject("eta_weights"))
 723   {
 724    fEtaWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("eta_weights"));
 725    if(!fEtaWeights){printf("\n WARNING (FQD): !fEtaWeights !!!!\n");exit(0);}
 726    if(TMath::Abs(fEtaWeights->GetBinWidth(1)-fEtaBinWidth)>pow(10.,-6.))
 727    {
 728     cout<<endl;
 729     cout<<"WARNING (FQD): Inconsistent binning in histograms for eta-weights throughout the code."<<endl;
 730     cout<<endl;
 731     //exit(0);
 732    }
 733   } else
 734     {
 735      cout<<"WARNING: fUseEtaWeights && fWeightsList->FindObject(\"eta_weights\") is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
 736      exit(0);
 737     }
 738  } // end of if(fUseEtaWeights)
 739
 740 } // end of AliFlowAnalysisWithFittingQDistribution::BookAndFillWeightsHistograms()
 741
 742 //================================================================================================================================
 743
 744 void AliFlowAnalysisWithFittingQDistribution::AccessConstants()
 745 {
 746  // access needed common constants from AliFlowCommonConstants
 747
 748  fnBinsPhi = AliFlowCommonConstants::GetMaster()->GetNbinsPhi();
 749  fPhiMin = AliFlowCommonConstants::GetMaster()->GetPhiMin();
 750  fPhiMax = AliFlowCommonConstants::GetMaster()->GetPhiMax();
 751  if(fnBinsPhi) fPhiBinWidth = (fPhiMax-fPhiMin)/fnBinsPhi;
 752  fnBinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 753  fPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
 754  fPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
 755  if(fnBinsPt) fPtBinWidth = (fPtMax-fPtMin)/fnBinsPt;
 756  fnBinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 757  fEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
 758  fEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
 759  if(fnBinsEta) fEtaBinWidth = (fEtaMax-fEtaMin)/fnBinsEta;
 760
 761 } // end of void AliFlowAnalysisWithFittingQDistribution::AccessConstants()
 762
 763 //================================================================================================================================
 764
 765 void AliFlowAnalysisWithFittingQDistribution::BookEverythingForDistributions()
 766 {
 767  // Book all objects relevant for fitting of q-distributions.
 768
 769  // Flags:
 770  TString sigmaFlag[2] = {"#sigma^{2} not fitted","#sigma^{2} fitted"};
 771  // q-distribution:
 772  TString fqDistributionName = "fqDistribution";
 773  fqDistributionName += fAnalysisLabel->Data();
 774  fqDistribution = new TH1D(Form("%s",fqDistributionName.Data()),"q-distribution",fqNbins,fqMin,fqMax);
 775  fqDistribution->SetXTitle(Form("q_{%d}=|Q_{%d}|/#sqrt{M}",fHarmonic,fHarmonic));
 776  fqDistribution->SetYTitle("Counts");
 777  fHistList->Add(fqDistribution);
 778  // q-distribution vs multiplicity:
 779  if(fStoreqDistributionVsMult)
 780  {
 781   TString fqDistributionVsMultName = "fqDistributionVsMult";
 782   fqDistributionVsMultName += fAnalysisLabel->Data();
 783   fqDistributionVsMult = new TH2D(Form("%s",fqDistributionVsMultName.Data()),"q-distribution vs M",fqNbins,fqMin,fqMax,fnBinsMult,fMinMult,fMaxMult);
 784   fqDistributionVsMult->GetXaxis()->SetTitle(Form("q_{%d}=|Q_{%d}|/#sqrt{M}",fHarmonic,fHarmonic));
 785   if(fMultiplicityIs==AliFlowCommonConstants::kRP)
 786   {
 787    fqDistributionVsMult->GetYaxis()->SetTitle("# RPs");
 788   } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
 789     {
 790      fqDistributionVsMult->GetYaxis()->SetTitle("Reference multiplicity (from ESD)");
 791     } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
 792       {
 793        fqDistributionVsMult->GetYaxis()->SetTitle("# POIs");
 794       }
 795   fqDistributionVsMult->GetZaxis()->SetTitle("Counts");
 796   fHistList->Add(fqDistributionVsMult);
 797  } // end of if(fStoreqDistributionVsMult)
 798  // Sum of particle weights:
 799  TString fSumOfParticleWeightsName = "fSumOfParticleWeights";
 800  fSumOfParticleWeightsName += fAnalysisLabel->Data();
 801  fSumOfParticleWeights = new TH1D(Form("%s",fSumOfParticleWeightsName.Data()),"Sum of particle weights",10000,0,100000); // (to be improved - harwired limits and number of bins)
 802  fSumOfParticleWeights->SetXTitle("#sum_{i=1}^{N} w_{i}");
 803  fSumOfParticleWeights->GetXaxis()->SetTitleSize(0.03);
 804  fSumOfParticleWeights->SetYTitle("Counts");
 805  fHistList->Add(fSumOfParticleWeights);
 806  // Final results for integrated flow:
 807  TString fIntFlowName = "fIntFlowFQD";
 808  fIntFlowName += fAnalysisLabel->Data();
 809  // sigma^2:
 810  TString fSigma2Name = "fSigma2";
 811  fSigma2Name += fAnalysisLabel->Data();
 812  // chi^2:
 813  TString fChi2Name = "fChi2";
 814  fChi2Name += fAnalysisLabel->Data();
 815  // fitting function:
 816  TString fittingFunctionName = "fFittingFunction";
 817  fittingFunctionName += fAnalysisLabel->Data();
 818
 819  for(Int_t f=0;f<2;f++) // sigma^2 not fitted (0) or fitted (1)
 820  {
 821   // final results for integrated flow:
 822   fIntFlow[f] = new TH1D(Form("%s, %s",fIntFlowName.Data(),sigmaFlag[f].Data()),Form("Reference Flow (%s)",sigmaFlag[f].Data()),1,0,1);
 823   fIntFlow[f]->SetLabelSize(0.08);
 824   fHistList->Add(fIntFlow[f]);
 825   // sigma^2:
 826   fSigma2[f] = new TH1D(Form("%s, %s",fSigma2Name.Data(),sigmaFlag[f].Data()),Form("#sigma^{2} (%s)",sigmaFlag[f].Data()),1,0,1);
 827   fSigma2[f]->SetLabelSize(0.08);
 828   (fSigma2[f]->GetXaxis())->SetBinLabel(1,"#sigma^{2}");
 829   fHistList->Add(fSigma2[f]);
 830   // chi^2:
 831   fChi2[f] = new TH1D(Form("%s, %s",fChi2Name.Data(),sigmaFlag[f].Data()),Form("#chi^{2} (%s)",sigmaFlag[f].Data()),1,0,1);
 832   fChi2[f]->SetLabelSize(0.08);
 833   (fChi2[f]->GetXaxis())->SetLabelOffset(0.01);
 834   (fChi2[f]->GetXaxis())->SetBinLabel(1,"#chi^{2}");
 835   fHistList->Add(fChi2[f]);
 836   // fitting functions:
 837   fFittingFunction[f] = new TF1(Form("%s, %s",fittingFunctionName.Data(),sigmaFlag[f].Data()),"[2]*(x/[1])*exp(-(x*x+[0]*[0])/(2.*[1]))*TMath::BesselI0(x*[0]/[1])");
 838   fHistList->Add(fFittingFunction[f]);
 839  } // end of for(Int_t f=0;f<2;f++) // sigma^2 not fitted or fitted
 840  // Book profile fFittingParameters which will hold all fitting parameters:
 841  TString fFittingParametersName = "fFittingParameters";
 842  fFittingParametersName += fAnalysisLabel->Data();
 843  fFittingParameters = new TProfile(fFittingParametersName.Data(),"Parameters for fitting q-distribution",12,0,12);
 844  fFittingParameters->SetLabelSize(0.05);
 845  fFittingParameters->GetXaxis()->SetBinLabel(1,"treshold");
 846  fFittingParameters->GetXaxis()->SetBinLabel(2,Form("starting v_{%d}",fHarmonic));
 847  fFittingParameters->GetXaxis()->SetBinLabel(3,Form("min. v_{%d}",fHarmonic));
 848  fFittingParameters->GetXaxis()->SetBinLabel(4,Form("max. v_{%d}",fHarmonic));
 849  fFittingParameters->GetXaxis()->SetBinLabel(5,"starting #sigma^{2}");
 850  fFittingParameters->GetXaxis()->SetBinLabel(6,"min. #sigma^{2}");
 851  fFittingParameters->GetXaxis()->SetBinLabel(7,"max. #sigma^{2}");
 852  fFittingParameters->GetXaxis()->SetBinLabel(8,"Final result is from #sigma^{2} fitted?");
 853  fFittingParameters->GetXaxis()->SetBinLabel(9,"Print results on the screen?");
 854  fFittingParameters->GetXaxis()->SetBinLabel(10,"fStoreqDistributionVsMult");
 855  fFittingParameters->GetXaxis()->SetBinLabel(11,"fMultiplicityIs");
 856  fFittingParameters->GetXaxis()->SetBinLabel(12,"fDoFit");
 857  fHistList->Add(fFittingParameters);
 858
 859 } // end of void AliFlowAnalysisWithFittingQDistribution::BookEverythingForDistributions()
 860
 861 //================================================================================================================================
 862
 863 void AliFlowAnalysisWithFittingQDistribution::DoFit(Bool_t sigma2Fitted)
 864 {
 865  // Do the final fit of q-distribution.
 866
 867  Int_t s2F = (Int_t)(sigma2Fitted); // shortcut
 868  Double_t AvM = fSumOfParticleWeights->GetMean(1); // average multiplicity
 869  //Int_t nEvts = (Int_t)fSumOfParticleWeights->GetEntries(); // number of events:
 870
 871  // Start fitting from the bin with at least fTreshold entries,
 872  // finish fitting at the bin with at least fTreshold entries:
 873  Int_t binMin = fqDistribution->FindFirstBinAbove(fTreshold);
 874  Int_t binMax = fqDistribution->FindLastBinAbove(fTreshold);
 875  Double_t binWidth = fqDistribution->GetBinWidth(4); // assuming that all bins have the same width
 876  if(TMath::Abs(binWidth) < 1.e-44)
 877  {
 878   cout<<endl;
 879   cout<<"WARNING (FQD): binWidth == 0 in AFAWFQD::DoFit()"<<endl;
 880   cout<<endl;
 881   exit(0);
 882  }
 883  Double_t qmin = (binMin-1)*binWidth;
 884  Double_t qmax = (binMax)*binWidth;
 885  Double_t ent = 0.; // number of entries between binMin and binMax:
 886  for(Int_t b=binMin;b<=binMax;b++)
 887  {
 888   ent += fqDistribution->GetBinContent(b);
 889  }
 890
 891  Double_t norm = binWidth*ent; // norm (assuming that all bins have the same width)
 892  // Fitting function:
 893  fFittingFunction[s2F]->SetRange(qmin,qmax);
 894  fFittingFunction[s2F]->SetParNames("v*sqrt{sum of particle weights}","sigma^2","norm");
 895  fFittingFunction[s2F]->SetParameters(fvStart*pow(AvM,0.5),fSigma2Start,norm);
 896  fFittingFunction[s2F]->SetParLimits(0,fvMin*pow(AvM,0.5),fvMax*pow(AvM,0.5));
 897  if(s2F == 0)
 898  {
 899   fFittingFunction[s2F]->FixParameter(1,0.5);
 900  } else
 901    {
 902     fFittingFunction[s2F]->SetParLimits(1,fSigma2Min,fSigma2Max);
 903    }
 904  fFittingFunction[s2F]->FixParameter(2,norm);
 905  // Fitting is done here:
 906  // Remark: do fit only if # of entries > 50 - this is only a pragmatics fix to avoid TMinuit crash (to be improved)
 907  if(ent > 50)
 908  {
 909   fqDistribution->Fit(fFittingFunction[s2F]->GetName(),"NQ","",qmin,qmax);
 910  }
 911
 912  // Final results:
 913  Double_t v = 0.; // reference flow
 914  Double_t vError = 0.; // error of reference flow
 915  Double_t sigma2 = 0.; // sigma^2
 916  Double_t sigma2Error = 0.; // error of sigma^2
 917  Double_t chi2 = 0; // chi^2 from Minuit
 918  // Reference flow:
 919  if(AvM)
 920  {
 921   v = fFittingFunction[s2F]->GetParameter(0)/pow(AvM,0.5);
 922   vError = fFittingFunction[s2F]->GetParError(0)/pow(AvM,0.5);
 923   fIntFlow[s2F]->SetBinContent(1,v); // s2F is shortcut for "sigma^2 fitted"
 924   fIntFlow[s2F]->SetBinError(1,vError); // s2F is shortcut for "sigma^2 fitted"
 925  } else
 926    {
 927     cout<<endl;
 928     cout<<"WARNING (FQD): AvM == 0 in AFAWFQD::DoFit()"<<endl;
 929     cout<<endl;
 930    }
 931  // sigma^2::
 932  if(s2F == 0) // sigma^2 not fitted, but fixed to 0.5
 933  {
 934   sigma2 = 0.5;
 935   fSigma2[0]->SetBinContent(1,sigma2);
 936   fSigma2[0]->SetBinError(1,0.);
 937  } else // sigma^2 fitted
 938    {
 939     sigma2 = fFittingFunction[s2F]->GetParameter(1);
 940     sigma2Error = fFittingFunction[s2F]->GetParError(1);
 941     fSigma2[1]->SetBinContent(1,sigma2);
 942     fSigma2[1]->SetBinError(1,sigma2Error);
 943    }
 944  // chi^2:
 945  chi2 = fFittingFunction[s2F]->GetChisquare();
 946  fChi2[s2F]->SetBinContent(1,chi2);
 947
 948 } // end of void AliFlowAnalysisWithFittingQDistribution::DoFit()
 949
 950 //================================================================================================================================
 951
 952 void AliFlowAnalysisWithFittingQDistribution::FillCommonHistResults(Bool_t sigma2Fitted)
 953 {
 954  // Fill common result histogram for reference flow and resolution.
 955
 956  // Remark: by default the result obtained with sigma^2 fitted is being stored.
 957  // In order to store the result obtained with sigma^2 fixed use a setter SetFinalResultIsFromSigma2Fitted(kFALSE).
 958
 959  Int_t s2F = (Int_t)sigma2Fitted; // shortcut
 960
 961  // Reference flow:
 962  Double_t v = fIntFlow[s2F]->GetBinContent(1);
 963  Double_t vError = fIntFlow[s2F]->GetBinError(1);
 964  fCommonHistsResults->FillIntegratedFlow(v,vError);
 965  // Resolution:
 966  Double_t AvM = fSumOfParticleWeights->GetMean(1);
 967  Double_t chi2 = AvM*pow(v,2.); // chi^2
 968  if(chi2>=0.)
 969  {
 970   fCommonHistsResults->FillChi(pow(chi2,0.5));
 971  }
 972
 973 } // end of void AliFlowAnalysisWithFittingQDistribution::FillCommonHistResultsIntFlow(Bool_t useParticleWeights, Bool_t sigma2NotFixed)
 974
 975 //================================================================================================================================
 976
 977 void AliFlowAnalysisWithFittingQDistribution::PrintOnTheScreen()
 978 {
 979  // Print the final results on the screen.
 980
 981  // shortcut for the harmonic:
 982  Int_t n = 0;
 983  if(fCommonHists->GetHarmonic())
 984  {
 985   n = (Int_t)(fCommonHists->GetHarmonic())->GetBinContent(1);
 986  }
 987
 988  // printing:
 989  cout<<endl;
 990  cout<<"***************************************"<<endl;
 991  cout<<"***************************************"<<endl;
 992  cout<<"      reference flow by fitting "<<endl;
 993  cout<<"           q-distribution:      "<<endl;
 994  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
 995  {
 996   cout<<"           (with weights)       "<<endl;
 997  } else
 998    {
 999     cout<<"          (without weights)       "<<endl;
1000    }
1001  cout<<endl;
1002  cout<<"1.) sigma^2 not fitted: "<<endl;
1003  cout<<"  v_"<<n<<"{FQD} = "<<fIntFlow[0]->GetBinContent(1)<<" +/- "<<fIntFlow[0]->GetBinError(1)<<endl;
1004  cout<<"  sigma^2 = 0.5 +/- 0 "<<endl;
1005  cout<<"  chi^2 = "<<fChi2[0]->GetBinContent(1)<<" (Minuit)"<<endl;
1006  cout<<endl;
1007  if(TMath::Abs(fIntFlow[0]->GetBinContent(1)-fvMin)<1.e-10 ||
1008     TMath::Abs(fIntFlow[0]->GetBinContent(1)-fvMax)<1.e-10)
1009  {
1010   cout<<"  WARNING: value of v_"<<n<<"{FQD}"<<" is on the boundary"<<endl;
1011   cout<<"           of fitting interval. Redo the fit."<< endl;
1012   cout<<endl;
1013  }
1014  cout<<"2.) sigma^2 fitted: "<<endl;
1015  cout<<"  v_"<<n<<"{FQD} = "<<fIntFlow[1]->GetBinContent(1)<<" +/- "<<fIntFlow[1]->GetBinError(1)<<endl;
1016  cout<<"  sigma^2 = "<<fSigma2[1]->GetBinContent(1)<<" +/- "<<fSigma2[1]->GetBinError(1)<<endl;
1017  cout<<"  chi^2 = "<<fChi2[1]->GetBinContent(1)<<" (Minuit)"<<endl;
1018  cout<<endl;
1019  if(TMath::Abs(fIntFlow[0]->GetBinContent(1)-fvMin)<1.e-10 ||
1020     TMath::Abs(fIntFlow[0]->GetBinContent(1)-fvMax)<1.e-10)
1021  {
1022   cout<<"  WARNING: value of v_"<<n<<"{FQD}"<<" is on the boundary"<<endl;
1023   cout<<"           of fitting interval. Redo the fit."<< endl;
1024   cout<<endl;
1025  }
1026  if(TMath::Abs(fSigma2[1]->GetBinContent(1)-fSigma2Min)<1.e-10 ||
1027     TMath::Abs(fSigma2[1]->GetBinContent(1)-fSigma2Max)<1.e-10)
1028  {
1029   cout<<"  WARNING: value of sigma^2 is on the boundary"<<endl;
1030   cout<<"           of fitting interval. Redo the fit."<< endl;
1031   cout<<endl;
1032  }
1033  cout<<"      nEvts = "<<fSumOfParticleWeights->GetEntries()<<", AvM = "<<fSumOfParticleWeights->GetMean()<<endl;
1034  cout<<endl;
1035  cout<<"***************************************"<<endl;
1036  cout<<"***************************************"<<endl;
1037  cout<<endl;
1038
1039 } // end of void AliFlowAnalysisWithFittingQDistribution::PrintOnTheScreen()
1040
1041 //================================================================================================================================
1042
1043 void AliFlowAnalysisWithFittingQDistribution::StoreFittingParameters()
1044 {
1045  // Store fitting parameters for the fit of q-distribution in profile fFittingParameters.
1046
1047  if(!fFittingParameters)
1048  {
1049   cout<<endl;
1050   cout<<"WARNING (FQD): fFittingParameters is NULL in AFAWFQD::SFP() !!!!"<<endl;
1051   cout<<endl;
1052   exit(0);
1053  }
1054  fFittingParameters->Reset();
1055  fFittingParameters->Fill(0.5,fTreshold);
1056  fFittingParameters->Fill(1.5,fvStart);
1057  fFittingParameters->Fill(2.5,fvMin);
1058  fFittingParameters->Fill(3.5,fvMax);
1059  fFittingParameters->Fill(4.5,fSigma2Start);
1060  fFittingParameters->Fill(5.5,fSigma2Min);
1061  fFittingParameters->Fill(6.5,fSigma2Max);
1062  fFittingParameters->Fill(7.5,fFinalResultIsFromSigma2Fitted);
1063  fFittingParameters->Fill(8.5,fPrintOnTheScreen);
1064  fFittingParameters->Fill(9.5,fStoreqDistributionVsMult);
1065  // which multiplicity was used:
1066  if(fMultiplicityIs==AliFlowCommonConstants::kRP) // # of Reference Particles
1067  {
1068   fFittingParameters->Fill(10.5,0); // 0 = # of Reference Particles
1069  } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
1070    {
1071     fFittingParameters->Fill(10.5,1); // 1 = ref. mult. from ESD
1072    } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
1073      {
1074       fFittingParameters->Fill(10.5,2); // 2 = # of Particles of Interest
1075      }
1076  fFittingParameters->Fill(11.5,fDoFit);
1077
1078 } // end of void AliFlowAnalysisWithFittingQDistribution::StoreFittingParameters()
1079
1080 //================================================================================================================================
1081
1082 void AliFlowAnalysisWithFittingQDistribution::AccessFittingParameters()
1083 {
1084  // Access fitting parameters for the fit of q-distribution.
1085
1086  fTreshold = fFittingParameters->GetBinContent(1);
1087  fvStart = fFittingParameters->GetBinContent(2);
1088  fvMin = fFittingParameters->GetBinContent(3);
1089  fvMax = fFittingParameters->GetBinContent(4);
1090  fSigma2Start = fFittingParameters->GetBinContent(5);
1091  fSigma2Min = fFittingParameters->GetBinContent(6);
1092  fSigma2Max = fFittingParameters->GetBinContent(7);
1093  fFinalResultIsFromSigma2Fitted = (Bool_t)fFittingParameters->GetBinContent(8);
1094  fPrintOnTheScreen = (Bool_t)fFittingParameters->GetBinContent(9);
1095  fStoreqDistributionVsMult = (Bool_t)fFittingParameters->GetBinContent(10);
1096  fDoFit = (Bool_t)fFittingParameters->GetBinContent(12);
1097
1098 } // end of void AliFlowAnalysisWithFittingQDistribution::AccessFittingParameters()