HBTAN/AliHBTCorrectCorrelFctn.cxx

   1 #include "AliHBTCorrectCorrelFctn.h"
   2 //____________________
   3 ///////////////////////////////////////////////////////
   4 //                                                   //
   5 // AliHBTCorrectQInvCorrelFctn                       //
   6 //                                                   //
   7 // Class for calculating Q Invariant correlation     //
   8 // taking to the account resolution of the           //
   9 // detector and coulomb effects.                     //
  10 // Implemented on basis of STAR procedure            //
  11 // implemented and described by Michael A. Lisa      //
  12 //                                                   //
  13 // Piotr.Skowronski@cern.ch                          //
  14 // http://alisoft.cern.ch/people/skowron/analyzer    //
  15 //                                                   //
  16 ///////////////////////////////////////////////////////
  17
  18 //Parameters fit from pi+ pi+ resolution analysis for pair with qinv<50MeV
  19 // chi2/NFD = 97/157
  20 // FCN=98.0971 FROM MIGRAD    STATUS=CONVERGED     332 CALLS         333 TOTAL
  21 //   EDM=2.13364e-12    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   2.4 per cent
  22 //  EXT PARAMETER                                   STEP         FIRST
  23 //  NO.   NAME         VALUE          ERROR         SIZE      DERIVATIVE
  24 //   1   fThetaA     2.72546e-03   1.43905e-04  -0.00000e+00   5.54375e-02
  25 //   2   fThetaB     1.87116e-04   5.11862e-05   0.00000e+00   3.66500e-01
  26 //   3  fThetaAlpha -2.36868e+00   1.83230e-01   0.00000e+00  -7.01301e-05
  27
  28 // FCN=120.603 FROM MIGRAD    STATUS=CONVERGED     117 CALLS         118 TOTAL
  29 //  EDM=2.5153e-12    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   2.4 per cent
  30 //  EXT PARAMETER                                   STEP         FIRST
  31 //  NO.   NAME         VALUE          ERROR         SIZE      DERIVATIVE
  32 //   1   fPhiA       1.93913e-03   1.31059e-04  -0.00000e+00  -5.87280e-02
  33 //   2   fPhiA       2.48687e-04   5.41251e-05   0.00000e+00  -1.87361e-01
  34 //   3  fPhiAlpha   -2.22649e+00   1.44503e-01   0.00000e+00   8.97538e-06
  35
  36 #include <TH1.h>
  37 #include <TH3.h>
  38 #include <TF1.h>
  39 #include <TRandom.h>
  40 ClassImp(AliHBTCorrectQInvCorrelFctn)
  41
  42 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title):
  43   AliHBTOnePairFctn1D(name,title),
  44   fMeasCorrelFctn(0x0),
  45   fMeasNumer(0x0),
  46   fMeasDenom(0x0),
  47   fSmearedNumer(0x0),
  48   fSmearedDenom(0x0),
  49   fDPtOverPtRMS(0.004),
  50   fThetaA(2.72e-03),
  51   fThetaB(1.87e-04),
  52   fThetaAlpha(-2.4),
  53   fPhiA(1.94e-03),
  54   fPhiB(2.5e-04),
  55   fPhiAlpha(-2.2),
  56   fR2(0.0),
  57   fLambda(0.0),
  58   fRConvergenceTreshold(0.3),
  59   fLambdaConvergenceTreshold(0.05)
  60 {
  61
  62 }
  63 /******************************************************************/
  64
  65 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title):
  66   AliHBTOnePairFctn1D(name,title,
  67                       measqinv->GetNbinsX(),
  68           measqinv->GetXaxis()->GetXmax(),
  69           measqinv->GetXaxis()->GetXmin()),
  70   fMeasCorrelFctn(measqinv),
  71   fMeasNumer(0x0),
  72   fMeasDenom(0x0),
  73   fSmearedNumer(0x0),
  74   fSmearedDenom(0x0),
  75   fDPtOverPtRMS(0.004),
  76   fThetaA(2.72e-03),
  77   fThetaB(1.87e-04),
  78   fThetaAlpha(-2.4),
  79   fPhiA(1.94e-03),
  80   fPhiB(2.5e-04),
  81   fPhiAlpha(-2.2),
  82   fR2(0.0),
  83   fLambda(0.0),
  84   fRConvergenceTreshold(0.3),
  85   fLambdaConvergenceTreshold(0.05)
  86 {
  87
  88 }
  89 /******************************************************************/
  90
  91 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title,
  92                     Int_t nbins, Float_t maxXval, Float_t minXval):
  93   AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval),
  94   fMeasCorrelFctn(0x0),
  95   fMeasNumer(0x0),
  96   fMeasDenom(0x0),
  97   fSmearedNumer(0x0),
  98   fSmearedDenom(0x0),
  99   fDPtOverPtRMS(0.004),
 100   fThetaA(2.72e-03),
 101   fThetaB(1.87e-04),
 102   fThetaAlpha(-2.4),
 103   fPhiA(1.94e-03),
 104   fPhiB(2.5e-04),
 105   fPhiAlpha(-2.2),
 106   fR2(0.0),
 107   fLambda(0.0),
 108   fRConvergenceTreshold(0.3),
 109   fLambdaConvergenceTreshold(0.05)
 110 {
 111 }
 112 /******************************************************************/
 113
 114 AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
 115 {
 116   delete fMeasCorrelFctn;
 117   delete fSmearedNumer;
 118   delete fSmearedDenom;
 119   delete fMeasNumer;
 120   delete fMeasDenom;
 121 }
 122 /******************************************************************/
 123
 124 void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
 125 {
 126   AliHBTFunction1D::BuildHistos(nbins,max,min);
 127
 128   TString numstr = fName + " Smeared Numerator";  //title and name of the numerator histogram
 129   TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram
 130
 131   fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
 132   fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
 133   fSmearedNumer->Sumw2();
 134   fSmearedDenom->Sumw2();
 135
 136   if (fMeasCorrelFctn == 0x0)
 137    {
 138      numstr = fName + " Measured Numerator";  //title and name of the numerator histogram
 139      denstr = fName + " Measured Denominator";//title and name of the denominator histogram
 140
 141      fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
 142      fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
 143      fMeasNumer->Sumw2();
 144      fMeasDenom->Sumw2();
 145    }
 146 }
 147 /******************************************************************/
 148
 149 void AliHBTCorrectQInvCorrelFctn::Init()
 150 {
 151   AliHBTOnePairFctn1D::Init();
 152   Info("Init","");
 153   fSmearedNumer->Reset();
 154   fSmearedDenom->Reset();
 155   if (fMeasNumer) fMeasNumer->Reset();
 156   if (fMeasDenom) fMeasDenom->Reset();
 157   fFittedR = -1.0;
 158   fFittedLambda = 0.0;
 159 }
 160 /******************************************************************/
 161
 162 void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
 163 {
 164   if (fMeasNumer == 0x0) return;
 165   pair = CheckPair(pair);
 166   if( pair == 0x0) return;
 167   fMeasNumer->Fill(pair->GetQInv());
 168 }
 169 /******************************************************************/
 170
 171 void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
 172 {
 173 //Process different events
 174   static AliHBTParticle part1, part2;
 175   static AliHBTPair smearedpair(&part1,&part2);
 176
 177   pair = CheckPair(pair);
 178   if( pair == 0x0) return;
 179   Double_t cc = GetCoulombCorrection(pair);
 180   Double_t qinv = pair->GetQInv();
 181
 182   //measured histogram -> if we are interested
 183   //only if fMeasCorrelFctn is not specified by user
 184   if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
 185
 186   Smear(pair,smearedpair);
 187   Double_t modelqinv = GetModelValue(qinv);
 188   //Ideal histogram
 189   fNumerator->Fill(qinv,modelqinv*cc);
 190   fDenominator->Fill(qinv,cc);
 191
 192   //Smeared histogram
 193   Double_t smearedqinv = smearedpair.GetQInv();
 194   fSmearedNumer->Fill(smearedqinv,modelqinv);
 195   Double_t smearedcc = GetCoulombCorrection(&smearedpair);
 196   fSmearedDenom->Fill(smearedqinv,smearedcc);
 197
 198 }
 199 /******************************************************************/
 200 void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
 201 {
 202 //Smears pair
 203   Smear(pair->Particle1(),smeared.Particle1());
 204   Smear(pair->Particle2(),smeared.Particle2());
 205   smeared.Changed();
 206 }
 207 /******************************************************************/
 208
 209 void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTParticle* part, AliHBTParticle* smeared)
 210 {
 211   Double_t sin2theta = TMath::Sin(part->Theta());
 212   sin2theta = sin2theta*sin2theta;
 213   Double_t pt = part->Pt();
 214
 215   double DpT_div_pT = gRandom->Gaus(0.0,fDPtOverPtRMS);
 216   double Dphi = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fPhiAlpha));
 217   double Dtheta = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fThetaAlpha));
 218
 219   Double_t smearedPx = part->Px()*(1.0+DpT_div_pT) - part->Py()*Dphi;
 220 //  fourmom.setX(px*(1.0+DpT_div_pT) - py*Dphi);
 221   Double_t smearedPy = part->Py()*(1.0+DpT_div_pT) - part->Px()*Dphi;
 222 //  fourmom.setY(py*(1.0+DpT_div_pT) + px*Dphi);
 223   Double_t smearedPz = part->Pz()*(1.0+DpT_div_pT) - pt*Dtheta/sin2theta;
 224 //  fourmom.setZ(pz*(1.0+DpT_div_pT) - pT*Dtheta/sin2theta);
 225
 226   Double_t mass2 = part->GetMass()*part->GetMass();
 227   Double_t E = mass2 + smearedPx*smearedPx +
 228                        smearedPy*smearedPy +
 229                        smearedPz*smearedPz;
 230
 231   smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(E));
 232 }
 233 /******************************************************************/
 234
 235 void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
 236 {
 237   fLambda = lambda;
 238   fR2 = r*r;
 239 }
 240 /******************************************************************/
 241 void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
 242 {
 243   //makes measured correlation function
 244   delete fMeasCorrelFctn;
 245   fMeasCorrelFctn = 0x0;
 246
 247   if (fMeasNumer&&fMeasDenom)
 248    {
 249      Double_t measscale = Scale(fMeasNumer,fMeasDenom);
 250      if (measscale == 0.0)
 251       {
 252         Error("GetResult","GetRatio for measured CF returned 0.0");
 253         return;
 254       }
 255      TString str = fName + "measured ratio";
 256      fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
 257      fMeasCorrelFctn->SetTitle(str.Data());
 258      fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
 259    }
 260 }
 261
 262 TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
 263 {
 264   //In case we don't have yet Measured Correlation Function
 265   //Try to get it
 266   //result is
 267   //        N[meas]   N[ideal]/D[ideal]
 268   //C(Q) =  ------- * -----------------
 269   //        D[meas]   N[smear]/D[smear]
 270
 271   TString str;
 272   if (fMeasCorrelFctn == 0x0) MakeMeasCF();
 273
 274   if (fMeasCorrelFctn == 0x0)
 275    {
 276      Error("GetResult",
 277            "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
 278      return 0x0;
 279    }
 280
 281   TH1D* ideal = (TH1D*)GetRatio(Scale());
 282   if (ideal == 0x0)
 283    {
 284      Error("GetResult","Ratio of ideal histograms is null");
 285      return 0x0;
 286    }
 287   str = fName + " smeared ratio";
 288   TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
 289   smearedCF->SetTitle(str.Data());
 290   Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
 291   smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
 292
 293   str = fName + " product meas ideal CF";
 294   TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
 295   measideal->Multiply(ideal,fMeasCorrelFctn);
 296
 297   str = fName + " Corrected Result";
 298   TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
 299   result->SetTitle(str.Data());
 300
 301   Double_t resultscale = Scale(measideal,smearedCF);
 302   result->Divide(measideal,smearedCF,resultscale);
 303   return result;
 304
 305 }
 306 /******************************************************************/
 307
 308 void AliHBTCorrectQInvCorrelFctn::Fit()
 309 {
 310 //fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))
 311 //where [0] is lambda
 312 //      [1] is radius
 313 //   0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
 314
 315   Info("Fit","Before fFittedLambda = %f",fFittedLambda);
 316   Info("Fit","Before fFittedR = %f",fFittedR);
 317   TH1D* result = (TH1D*)GetResult();
 318   if (result == 0x0)
 319    {
 320      Error("Fit","Can not get result");
 321      return;
 322    }
 323   TF1* fitfctn = new TF1("fitfctn","1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))");
 324   fitfctn->SetParameter(0,1.0);
 325   fitfctn->SetParameter(1,6.0);
 326   Float_t max = result->GetXaxis()->GetXmax();
 327   Info("Fit","Max is %f",max);
 328   result->Fit(fitfctn,"0","",0.008,max);
 329   fFittedLambda = fitfctn->GetParameter(0);
 330   fFittedR = fitfctn->GetParameter(1);
 331   Info("Fit","After fFittedLambda = %f",fFittedLambda);
 332   Info("Fit","After fFittedR = %f",fFittedR);
 333   delete fitfctn;
 334   delete result;
 335 }
 336 /******************************************************************/
 337
 338 Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
 339 {
 340   //check if fitting was performed
 341   if (fFittedR <= 0.0)
 342    {
 343      Fit();//if not do fit first
 344      if (fFittedR <= 0.0)
 345       {
 346         Error("IsConverged","Fitting failed");
 347         return kFALSE;
 348       }
 349    }
 350
 351   Double_t guessedR = TMath::Sqrt(fR2);
 352
 353   Info("IsConverged","Fitted   lambda            : %8f Fitted  Radius             : %8f",fFittedLambda,fFittedR);
 354   Info("IsConverged","Guessed  lambda            : %8f Guessed Radius             : %8f",fLambda,guessedR);
 355   Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f",
 356        fLambdaConvergenceTreshold,fRConvergenceTreshold);
 357
 358   if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) &&
 359        (TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
 360     {
 361       Info("IsConverged","Cnvergence reached");
 362       return kTRUE;
 363     }
 364   else
 365    {
 366       Info("IsConverged","Cnvergence NOT reached");
 367       return kFALSE;
 368    }
 369 }
 370 /******************************************************************/
 371
 372 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
 373 {
 374   if (fFittedR <= 0.0) Fit();
 375   return fFittedR;
 376 }
 377 /******************************************************************/
 378
 379 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
 380 {
 381   if (fFittedR <= 0.0) Fit();
 382   return fFittedLambda;
 383 }
 384 /******************************************************************/
 385
 386 void AliHBTCorrectQInvCorrelFctn::WriteAll()
 387 {
 388   Write();
 389   if (fMeasCorrelFctn) fMeasCorrelFctn->Write();
 390   if (fMeasNumer ) fMeasNumer->Write();
 391   if (fMeasDenom ) fMeasDenom->Write();
 392   if (fSmearedNumer) fSmearedNumer->Write();
 393   if (fSmearedDenom) fSmearedDenom->Write();
 394   if (fSmearedNumer && fSmearedDenom)
 395    {
 396     TString str = fName + " smeared ratio";
 397     TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
 398     smearedCF->SetTitle(str.Data());
 399     Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
 400     smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
 401    }
 402 }
 403
 404 /******************************************************************/
 405 /******************************************************************/
 406 /******************************************************************/
 407
 408 //____________________
 409 ///////////////////////////////////////////////////////
 410 //                                                   //
 411 // AliHBTCorrectQ3DCorrelFctn                        //
 412 //                                                   //
 413 // Class for calculating Q Invariant correlation     //
 414 // taking to the account resolution of the           //
 415 // detector and coulomb effects.                     //
 416 //                                                   //
 417 ///////////////////////////////////////////////////////
 418
 419
 420 AliHBTCorrectQ3DCorrelFctn::AliHBTCorrectQ3DCorrelFctn(const char* name, const char* title):
 421  AliHBTOnePairFctn3D(name,title),
 422  fMeasCorrelFctn(0x0),
 423  fSmearedNumer(0x0),
 424  fSmearedDenom(0x0),
 425  fMeasNumer(0x0),
 426  fMeasDenom(0x0)
 427 {
 428
 429 }
 430 /******************************************************************/
 431
 432 AliHBTCorrectQ3DCorrelFctn::~AliHBTCorrectQ3DCorrelFctn()
 433 {
 434  delete fMeasCorrelFctn;
 435  delete fSmearedNumer;
 436  delete fSmearedDenom;
 437  delete fMeasNumer;
 438  delete fMeasDenom;
 439 }
 440
 441 /******************************************************************/