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