1 /**************************************************************************
2 * Author: Andrey Ivanov. *
3 * Contributors are mentioned in the code where appropriate. *
5 * Permission to use, copy, modify and distribute this software and its *
6 * documentation strictly for non-commercial purposes is hereby granted *
7 * without fee, provided that the above copyright notice appears in all *
8 * copies and that both the copyright notice and this permission notice *
9 * appear in the supporting documentation. The authors make no claims *
10 * about the suitability of this software for any purpose. It is *
11 * provided "as is" without express or implied warranty. *
12 **************************************************************************/
13 //-------------------------------------------------------------------------
15 // This class include into LRC library for Long-Range Correlation analysis
16 // it is base class for NN, PtN, PtPt
17 // implements base methods for thees classes
18 // Origin: Petr Naumenko, SPbSU-CERN, Petr.Naoumenko@cern.ch ,
19 // Andrey Ivanov (SPbSU-CERN), Igor Altsebeev (SPbSU-CERN)
20 //-------------------------------------------------------------------------
22 #include "AliLRCAnalysis.h"
24 #include "AliLRCFit.h"
28 #include "TPaveText.h"
35 ClassImp(AliLRCAnalysis)
37 /******************************************************
38 * AliLRCAnalysis class
39 ******************************************************/
41 AliLRCAnalysis::AliLRCAnalysis():TObject(), fPrAbs(new TH1D()), fPrRel(new TH1D()), fPrf(new TH1D()), fPrb(new TH1D()), fileHist(new TFile()), fdptb(.0), fEntries(0), fSx((char*)" "), fSy((char*)" "), fxFitMin(.0), fxFitMax(.0), fNsigma(2.), farel(.0), fbrel(.0), farelError(.0), fbrelError(.0), fXi2rel(.0), faabs(.0), fbabs(.0), faabsError(.0), fbabsError(.0), fXi2abs(.0),fFitMethod(0){
45 AliLRCAnalysis::AliLRCAnalysis(const AliLRCAnalysis& a):TObject(a),fPrAbs(a.fPrAbs), fPrRel(a.fPrRel), fPrf(a.fPrf), fPrb(a.fPrb), fileHist(a.fileHist), fdptb(a.fdptb), fEntries(a.fEntries), fSx(a.fSx), fSy(a.fSy), fxFitMin(a.fxFitMin), fxFitMax(a.fxFitMax), fNsigma(a.fNsigma), farel(a.farel), fbrel(a.fbrel), farelError(a.farelError), fbrelError(a.fbrelError), fXi2rel(a.fXi2rel), faabs(a.faabs), fbabs(a.fbabs),faabsError(a.faabsError), fbabsError(a.fbabsError), fXi2abs(a.fXi2abs),fFitMethod(a.fFitMethod){
49 AliLRCAnalysis& AliLRCAnalysis::operator= (const AliLRCAnalysis& a){
52 TObject::operator= (a);
57 fileHist = a.fileHist;
61 fEntries = a.fEntries;
62 fxFitMin = a.fxFitMin;
63 fxFitMax = a.fxFitMax;
67 farelError = a.farelError;
68 fbrelError = a.fbrelError;
72 faabsError = a.faabsError;
73 fbabsError = a.fbabsError;
75 fFitMethod = a.fFitMethod;
80 AliLRCAnalysis::~AliLRCAnalysis() {
91 void AliLRCAnalysis::CreateHist(char *name, char *nameAbs, char *nameRel, char *atitleF, char *atitleB, char *rtitleF, char *rtitleB, TH2D* sourceHist) {
92 //Create absolute and relation var histogramm
93 TProfile* profX = (TProfile*) sourceHist->ProfileX(name, 1, sourceHist->GetNbinsY());
94 fEntries = (int) sourceHist->GetEntries();
95 fPrAbs = new TH1D(nameAbs, profX->GetTitle(), profX->GetXaxis()->GetNbins(), profX->GetXaxis()->GetXmin(), profX->GetXaxis()->GetXmax());
96 fPrAbs->SetOption("E");
97 fPrRel = new TH1D(nameRel, profX->GetTitle(), profX->GetXaxis()->GetNbins(), profX->GetXaxis()->GetXmin()/sourceHist->ProjectionX()->GetMean(), profX->GetXaxis()->GetXmax()/sourceHist->ProjectionX()->GetMean());
98 fPrRel->SetOption("E");
99 fPrAbs->GetXaxis()->SetTitle(atitleF);
100 fPrAbs->GetYaxis()->SetTitle(atitleB);
101 fPrRel->GetXaxis()->SetTitle(rtitleF);
102 fPrRel->GetYaxis()->SetTitle(rtitleB);
103 fPrf = (TH1D*) sourceHist->ProjectionX();
104 fPrb = (TH1D*) sourceHist->ProjectionY();
105 fPrf->GetXaxis()->SetTitle(atitleF);
106 fPrf->GetYaxis()->SetTitle("Tracks");
107 fPrb->GetXaxis()->SetTitle(atitleB);
108 fPrb->GetYaxis()->SetTitle("Tracks");
111 double mnf = fPrf->GetMean();
112 double rms = fPrf->GetRMS();
113 fxFitMin = mnf-fNsigma*rms;//sqrt(rms);
114 fxFitMax = mnf+fNsigma*rms;//sqrt(rms);//mnf+2*sqrt(rms);
119 void AliLRCAnalysis::SetBinsRange(int binMin, int binMax){
122 int n=h->GetNbinsX();
123 fxFitMin = h->GetXaxis()->GetXmin();
124 fxFitMax = h->GetXaxis()->GetXmax();
125 double df = (fxFitMax-fxFitMin)/n;
126 for(int i=1; i<=n; i++)
127 if(h->GetBinContent(i) != 0){
128 fxFitMin += (i + binMin) * df;
131 for(int i=1; i<=n; i++)
132 if(h->GetBinContent(n-i) != 0){
133 fxFitMax -= (i + binMax) * df;
138 bool AliLRCAnalysis::SetFitRange(double xMin, double xMax){
143 this->fxFitMin = xMin;
144 this->fxFitMax = xMax;
147 bool AliLRCAnalysis::SetFitRangeMin(double xMin){
148 //Set the fit range min
153 bool AliLRCAnalysis::SetFitRangeMax(double xMax){
154 //Set the fit range min
155 //this->fxFitMax = xMax;
158 double maxBorder = 0.;
160 for ( int binI = h->GetNbinsX(); binI > 0; binI-- )
162 if ( h->GetBinContent(binI) != 0 )
164 maxBorder = h->GetBinLowEdge(binI+1) ;
168 fxFitMax = TMath::Max( fxFitMin, TMath::Min( xMax, maxBorder ) );
173 void AliLRCAnalysis::SetFullFitRange(){
174 //Set fitting on full range
176 fxFitMin = h->GetXaxis()->GetXmin();
178 for ( int binI = h->GetNbinsX(); binI > 0; binI-- )
180 if ( h->GetBinContent(binI) != 0 )
182 fxFitMax = h->GetBinLowEdge(binI+1) ;
189 void AliLRCAnalysis::SetXmin(double xMin){
193 void AliLRCAnalysis::SetXmax(double xMax){
196 void AliLRCAnalysis::SetNsigma(double nSigma){
197 //Sets fiting range in forward value RMSs
199 double mnf = fPrf->GetMean();
200 double rms = fPrf->GetRMS();
201 fxFitMin = mnf-fNsigma*rms;
202 fxFitMax = mnf+fNsigma*rms;
204 // cout << "### Fit params: N sigma is " << fNsigma
205 // << ", mean = " << mnf << ", rms = " << rms << ", minFitX = " << fxFitMin << ", maxFitX = " << fxFitMax << endl;
208 void AliLRCAnalysis::SetFitMethod(int id){
212 double AliLRCAnalysis::GetFitXmin() const
217 double AliLRCAnalysis::GetFitXmax() const
221 double AliLRCAnalysis::GetArel() const {
225 double AliLRCAnalysis::GetBrel() const {
229 double AliLRCAnalysis::GetArelError() const {
233 double AliLRCAnalysis::GetBrelError() const {
237 double AliLRCAnalysis::GetXi2rel() const {
241 double AliLRCAnalysis::GetAabs() const {
245 double AliLRCAnalysis::GetBabs() const {
249 double AliLRCAnalysis::GetAabsError() const {
253 double AliLRCAnalysis::GetBabsError() const {
256 double AliLRCAnalysis::GetXi2abs() const {
260 void AliLRCAnalysis::Calculate(){
262 double mnf = fPrf->GetMean();
265 AliLRCFit *fit1 = new AliLRCFit(fPrRel, fxFitMin/mnf, fxFitMax/mnf, 1.);
266 AliLRCFit *fit2 = new AliLRCFit(fPrAbs, fxFitMin, fxFitMax, 0.);
267 farel = fit1->Geta();
268 fbrel = fit1->Getb();
270 faabs = fit2->Geta();
271 fbabs = fit2->Getb();
275 faabsError = fit2->Getda();
276 fbabsError = fit2->Getdb();
278 farelError = fit1->Getda();
279 fbrelError = fit1->Getdb();
283 faabsError = fit2->Getda1();
284 fbabsError = fit2->Getdb1();
286 farelError = fit1->Getda1();
287 fbrelError = fit1->Getdb1();
290 fXi2rel = fit1->Gethi2();
291 fXi2abs = fit2->Gethi2();
294 void AliLRCAnalysis::DrawAbs() {
295 //Draw abs var hist with ALL info
296 int * mas = new int [fgkPlotFlags];
297 for ( int i = 0; i < fgkPlotFlags; i++ )
299 DrawAbsPure( mas, 1 );
303 void AliLRCAnalysis::DrawAbs( const int * const mas ) {
304 //Draw abs var hist with REQUESTED BY ARRAY info
305 DrawAbsPure( mas, 1 );
308 void AliLRCAnalysis::DrawAbsPure( const int * const mDrawArray, bool drawPaveLabel ) {
310 // Draw abs var histrogram
311 DrawHist( mDrawArray, drawPaveLabel, faabs, fbabs, faabsError, fbabsError, fPrAbs, 0 );
314 void AliLRCAnalysis::DrawRel() {
315 //Draw rel var hist with ALL info
316 int * mas = new int [fgkPlotFlags];
317 for ( int i = 0; i < fgkPlotFlags; i++ )
319 DrawRelPure( mas, 1 );
323 void AliLRCAnalysis::DrawRel( const int * const mas ) {
324 //Draw rel var hist with REQUESTED BY ARRAY info
325 DrawRelPure( mas, 1 );
328 void AliLRCAnalysis::DrawRelPure( const int * const mDrawArray, bool drawPaveLabel ) {
330 // Draw rel var histogram
331 DrawHist( mDrawArray, drawPaveLabel, farel, fbrel, farelError, fbrelError, fPrRel, 1 );
333 void AliLRCAnalysis::DrawHist( const int * const mDrawArray, bool drawPaveLabel, double aCoef, double bCoef,
334 double aCoefError, double bCoefError, TH1D* profToDraw, int histType ) {
335 // Method called by DrawRelPure or DrawAbsPure to draw corresponding LRC histo
336 double mnf = fPrf->GetMean();
337 double y1, y2, x1, x2;
340 //mnf = fPrf->GetMean();
341 profToDraw->SetStats(0);
342 //profToDraw->Fit("pol1", "", "", fxFitMin, fxFitMax);
343 //profToDraw->Fit("pol1");
344 //TF1 *fit1 = profToDraw->GetFunction("pol1");
348 if ( histType == 0 ) {
349 f1 = new TF1( "f1", "[0]+[1]*x", fxFitMin, fxFitMax);
350 f1->SetLineColor(kRed);
351 f1->SetParameters(aCoef,bCoef);
353 else if ( histType == 1 ) {
354 f1 = new TF1( "f1", "[0]+[1]*(x-1)", fxFitMin/mnf, fxFitMax/mnf);
355 f1->SetLineColor(kGreen);
356 f1->SetParameters(aCoef,bCoef);
358 //cout << " set draw params: a=" << aCoef << " b=" << bCoef << endl;
360 y1=profToDraw->GetBinContent(1)-profToDraw->GetBinError(1);
361 y2=profToDraw->GetBinContent(1)+profToDraw->GetBinError(1);
362 n=profToDraw->GetNbinsX();
366 if(profToDraw->GetBinContent(i)-profToDraw->GetBinError(i)<y1)
367 y1=profToDraw->GetBinContent(i)-profToDraw->GetBinError(i);
368 if(profToDraw->GetBinContent(i)+profToDraw->GetBinError(i)>y2)
369 y2=profToDraw->GetBinContent(i)+profToDraw->GetBinError(i);
372 profToDraw->DrawCopy();
373 if(f1) f1->DrawCopy("same");
376 x1 = profToDraw->GetXaxis()->GetXmin();
377 x2 = profToDraw->GetXaxis()->GetXmax();
382 for ( int j = 0; j < 9; j++)
383 if ( mDrawArray[j] ) nDatas++;
384 double aXshift = (x2-x1)/7;
385 double aYshift = (y2-y1)/20;
387 TPaveText *pt1 = new TPaveText(x1+(x2-x1)/2 + aXshift, y1+aYshift, x2-(x2-x1)/6 + aXshift, y1+(y2-y1)/3*2/9*nDatas + aYshift);
388 snprintf(str,50, "Events = %i", fEntries);
389 //sprintf(str, "Events = %i", fEntries);
390 if ( mDrawArray[0] ) pt1->AddText(str);
391 snprintf(str, 50,"a = %g #pm %g", GetRoundWithError( aCoef, aCoefError ), GetRoundWithPrecision(aCoefError, 2)); //fit1->GetParameter(0), fit1->GetParError(0));
392 //sprintf(str, "a = %g #pm %g", GetRoundWithError( aCoef, aCoefError ), GetRoundWithPrecision(aCoefError, 2)); //fit1->GetParameter(0), fit1->GetParError(0));
393 if ( mDrawArray[1] ) pt1->AddText(str);
394 snprintf(str, 50,"b = %g #pm %g", GetRoundWithError( bCoef, bCoefError ), GetRoundWithPrecision(bCoefError, 2)); //fit1->GetParameter(1), fit1->GetParError(1));
395 //sprintf(str, "b = %g #pm %g", GetRoundWithError( bCoef, bCoefError ), GetRoundWithPrecision(bCoefError, 2)); //fit1->GetParameter(1), fit1->GetParError(1));
396 if ( mDrawArray[2] ) pt1->AddText(str);
397 snprintf(str, 50,"#hat{#chi}^{2} = #chi^{2}/(n-2) = %g", GetRoundWithPrecision(fXi2abs, 3));
398 //sprintf(str, "#hat{#chi}^{2} = #chi^{2}/(n-2) = %g", GetRoundWithPrecision(fXi2abs, 3));
399 if ( mDrawArray[3] ) pt1->AddText(str);
400 snprintf(str,50, "<%s> = %g " , fSx, GetRoundWithPrecision(fPrf->GetMean(), 3));
401 //sprintf(str, "<%s> = %g " , fSx, GetRoundWithPrecision(fPrf->GetMean(), 3));
402 if ( mDrawArray[4] ) pt1->AddText(str);
404 snprintf(str,50, "<%s> = %g", fSy, GetRoundWithPrecision(fPrb->GetMean(),3));
405 //sprintf(str, "<%s> = %g", fSy, GetRoundWithPrecision(fPrb->GetMean(),3));
406 if ( mDrawArray[5] ) pt1->AddText(str);
408 snprintf(str, 50,"<<%s>> = %g " , fSx, GetRoundWithPrecision(fPrf->GetRMS(), 3));
409 //sprintf(str, "<<%s>> = %g " , fSx, GetRoundWithPrecision(fPrf->GetRMS(), 3));
410 if ( mDrawArray[6] ) pt1->AddText(str);
411 snprintf(str, 50,"<<%s>> = %g", fSy, GetRoundWithPrecision(fPrb->GetRMS(), 3));
412 //sprintf(str, "<<%s>> = %g", fSy, GetRoundWithPrecision(fPrb->GetRMS(), 3));
413 if ( mDrawArray[7] ) pt1->AddText(str);
416 snprintf(str,50, "d%s = %g", fSy, GetRoundWithPrecision(fdptb, 3));
417 //sprintf(str, "d%s = %g", fSy, GetRoundWithPrecision(fdptb, 3));
418 if ( mDrawArray[8] ) pt1->AddText(str);
421 pt1->SetTextAlign(12);
422 pt1->SetTextFont(42);
423 pt1->SetFillColor(4000);
424 //pt1->SetFillStyle(4100);
425 pt1->SetShadowColor(4000);
426 pt1->SetBorderSize(0);
428 pt1->DrawClone("same");
432 void AliLRCAnalysis::SetGraphics() const {
433 // Set root graph style
435 gStyle->SetCanvasColor(10);
436 gStyle->SetFrameFillColor(10);
437 gStyle->SetStatColor(0);
438 gStyle->SetPadColor(0);
441 double AliLRCAnalysis::Integral(TH2D* const source, Int_t nbin) const {
442 // calculate the integrall for x bin and y bins of 2d histogramm
444 for(Int_t i = 1; i<=source->GetNbinsY(); i++) {
445 sum += source->GetYaxis()->GetBinCenter(i)*source->GetBinContent(nbin, i);
450 void AliLRCAnalysis::SetErrors(TH2D* source, const char *name){
451 // Calculate errors for NN
452 TProfile* profX = (TProfile*) source->ProfileX(name, 1, source->GetNbinsY());
453 for(int i = 0; i < profX->GetNbinsX(); i++)
455 fPrAbs->SetBinContent(i, profX->GetBinContent(i));
456 if(fPrf->GetBinContent(i)!=0)
458 fPrAbs->SetBinError(i,sqrt(profX->GetBinContent(i)/fPrf->GetBinContent(i)));
460 fPrRel->SetBinContent(i, fPrAbs->GetBinContent(i)/fPrb->GetMean());
461 fPrRel->SetBinError(i,fPrAbs->GetBinError(i)/fPrb->GetMean());
465 void AliLRCAnalysis::SetErrors(TH2D* source, const char *name, double ptd, TH2D* nb){
466 //Calculate arrors for ptn and ptpt
467 TProfile* profX = (TProfile*) source->ProfileX(name, 1, source->GetNbinsY());
470 for(int i = 0; i < profX->GetNbinsX(); i++)
472 fPrAbs->SetBinContent(i, profX->GetBinContent(i));
473 if(fPrf->GetBinContent(i)!=0)
475 pt = profX->GetBinContent(i);
476 fPrAbs->SetBinError(i,ptd*sqrt(Integral(nb,i))/fPrf->GetBinContent(i));
478 fPrRel->SetBinContent(i, fPrAbs->GetBinContent(i)/fPrb->GetMean());
479 fPrRel->SetBinError(i,fPrAbs->GetBinError(i)/fPrb->GetMean());
484 void AliLRCAnalysis::SetErrors(TH2D* source, const char *name, double ptd, const TProfile* nb){
485 //Calculate arrors for ptn and ptpt
486 TProfile* profX = (TProfile*) source->ProfileX(name, 1, source->GetNbinsY());
489 for(int i = 0; i < profX->GetNbinsX(); i++)
491 fPrAbs->SetBinContent(i, profX->GetBinContent(i));
492 if(fPrf->GetBinContent(i)!=0)
494 pt = profX->GetBinContent(i);
495 fPrAbs->SetBinError(i,ptd*sqrt(nb->GetBinContent(i)/fPrf->GetBinContent(i)));
497 fPrRel->SetBinContent(i, fPrAbs->GetBinContent(i)/fPrb->GetMean());
498 fPrRel->SetBinError(i,fPrAbs->GetBinError(i)/fPrb->GetMean());
503 double AliLRCAnalysis::GetRoundWithError( double value, double error ) const {
504 //Rounding error and value with DEFAULT precision
505 return GetRoundValueErrorPrecision( value, error, 2 );
507 double AliLRCAnalysis::GetRoundWithError( double value, double error, int pres ) const {
508 //Rounding error and value with REQUESTED precision
509 return GetRoundValueErrorPrecision( value, error, pres );
511 double AliLRCAnalysis::GetRoundWithPrecision( double value, int pres ) const {
512 //Rounding error and value with requested precision
513 return GetRoundValueErrorPrecision( value, 0, pres );
516 double AliLRCAnalysis::GetRoundValueErrorPrecision( double value, double error, int pres ) const {
517 //Rounding error and value with requested precision
518 //value == value, error == error
519 //if single argument(value=0) - calculate without errors:
523 bool noError = false;
526 //cout << "Before rounding: " << value << " " << error << endl;
533 while ( ((int)error)%10 == 0 && i < 6 )
540 for ( int j = 0; j < pres; j++ )
547 if ( ((int)(error*10))%10 > 4 && ((int)(error*10))%10 != 9 ) //trouble: if we round 19 to 20 - zero disappeares!
549 error = (double)((int)error + adding)/order;
553 //cout << "After rounding: " << error << endl;
558 for ( int j = 0; j < i; j++ )
562 if ( ((int)(value*10))%10 > 4 && ((int)(value*10))%10 != 9 )
564 value = (double)((int)value + adding)/order;
566 //cout << "After rounding: " << value << " " << error << endl;
567 return value; //taking into account ERROR
571 TH1D* AliLRCAnalysis::GetAbsHisto() const
573 //Returns final histo in absolute variables
576 TH1D* AliLRCAnalysis::GetRelHisto() const
578 //Returns final histo in relative variables
582 TH1D* AliLRCAnalysis::GetForwardValueDist() const
584 //Returns destribution of value used in forward window
587 TH1D* AliLRCAnalysis::GetBackwardValueDist() const
589 //Returns destribution of value used in backward window