HBTAN/AliHBTCorrectQInvCorrelFctn.cxx

   1 #include "AliHBTCorrectQInvCorrelFctn.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://aliweb.cern.ch/people/skowron/analyzer    //
  15 //                                                   //
  16 ///////////////////////////////////////////////////////
  17
  18 // dPt/Pt
  19 // root [19] rms->Fit(f,"","",0.1,2)
  20 //  FCN=7.0017 FROM MIGRAD    STATUS=CONVERGED     126 CALLS         127 TOTAL
  21 //                      EDM=2.28804e-15    STRATEGY= 1      ERROR MATRIX ACCURATE
  22 //    EXT PARAMETER                                   STEP         FIRST
  23 //  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
  24 //    1  p0           5.78220e-03   3.14576e-05   4.97822e-09  -1.90059e-05
  25 //    2  p1           3.98063e-05   1.61877e-06   1.04380e-10   1.91454e-04
  26 //    3  p2          -2.78008e+00   2.13581e-02   1.66031e-06   3.16574e-06
  27 //    4  p3           5.07594e-04   4.79619e-05   1.29793e-08  -2.29242e-05
  28
  29
  30 // Phi
  31 // root [17] rms->Fit(f,"","",0.15,2.5)
  32 // Warning in <TH1D::Fit>: Abnormal termination of minimization.
  33 // FCN=33.4898 FROM MIGRAD    STATUS=FAILED         91 CALLS          92 TOTAL
  34 //                     EDM=1.19154e-15    STRATEGY= 1      ERR MATRIX APPROXIMATE
  35 //    EXT PARAMETER                APPROXIMATE        STEP         FIRST
  36 //   NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
  37 //    1  p0           5.87693e-04   5.04254e-06   2.49187e-09   5.84546e-04
  38 //    2  p1           2.16488e-06   3.68880e-07   6.41507e-11  -7.36564e-02
  39 //    3  p2          -3.10218e+00   1.01695e-01   2.00177e-05   7.54285e-07
  40 //    4  p3          -1.79892e-05   5.44067e-06   2.15870e-09   4.11441e-04
  41
  42
  43
  44 // Theta
  45 // root [14] rms->Fit(f,"","",0.1,3)
  46 //  FCN=8.9981 FROM MIGRAD    STATUS=CONVERGED      79 CALLS          80 TOTAL
  47 //                      EDM=3.03049e-17    STRATEGY= 1      ERR MATRIX NOT POS-DEF
  48 //  EXT PARAMETER                APPROXIMATE        STEP         FIRST
  49 //  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
  50 //     1  p0          -1.68773e-02   2.67644e-05   8.04770e-09   4.48079e-05
  51 //     2  p1           1.78440e-02   2.65467e-05   8.50867e-09   6.43012e-05
  52 //     3  p2          -5.26559e-02   5.06308e-04   2.28595e-07  -1.63963e-05
  53 //     4  p3           2.00940e-04   1.14440e-05   3.98737e-09   1.78198e-05
  54
  55
  56 #include <TH1.h>
  57 #include <TH3.h>
  58 #include <TF1.h>
  59 #include <TRandom.h>
  60 #include <AliAODParticle.h>
  61
  62
  63 ClassImp(AliHBTCorrectedCorrelFctn)
  64
  65 AliHBTCorrectedCorrelFctn::AliHBTCorrectedCorrelFctn():
  66   fDPtOverPtA(5.78220e-03),
  67   fDPtOverPtB(3.98063e-05),
  68   fDPtOverPtAlpha(-2.78008),
  69   fDPtOverPtC(5.07594e-04),
  70   fThetaA(5.87693e-04),
  71   fThetaB(2.16488e-06),
  72   fThetaAlpha(-3.10218e+00),
  73   fThetaC(-1.79892e-05),
  74   fPhiA(-1.68773e-02),
  75   fPhiB(1.78440e-02),
  76   fPhiAlpha(-5.26559e-02),
  77   fPhiC(2.00940e-04)
  78 {
  79   //ctor
  80 }
  81
  82 /******************************************************************/
  83 void AliHBTCorrectedCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
  84 {
  85 //Smears pair
  86   Smear(pair->Particle1(),smeared.Particle1());
  87   Smear(pair->Particle2(),smeared.Particle2());
  88   smeared.Changed();
  89 }
  90 /******************************************************************/
  91
  92 void AliHBTCorrectedCorrelFctn::Smear(AliVAODParticle* part, AliVAODParticle* smeared)
  93 {
  94  //Smears momenta
  95   Double_t sin2theta = TMath::Sin(part->Theta());
  96   sin2theta = sin2theta*sin2theta;
  97   Double_t pt = part->Pt();
  98
  99   Double_t sigmapt = fDPtOverPtA + fDPtOverPtB*TMath::Power(pt,fDPtOverPtAlpha) + fDPtOverPtC*pt;
 100   Double_t dPtDivPt = gRandom->Gaus(0.0,sigmapt);
 101   Double_t dphi = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fPhiAlpha) + fPhiC*pt);
 102   Double_t dtheta = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fThetaAlpha) +fThetaC*pt);
 103
 104   Double_t smearedPx = part->Px()*(1.0+dPtDivPt) - part->Py()*dphi;
 105 //  fourmom.setX(px*(1.0+dPtDivPt) - py*dphi);
 106   Double_t smearedPy = part->Py()*(1.0+dPtDivPt) - part->Px()*dphi;
 107 //  fourmom.setY(py*(1.0+dPtDivPt) + px*dphi);
 108   Double_t smearedPz = part->Pz()*(1.0+dPtDivPt) - pt*dtheta/sin2theta;
 109 //  fourmom.setZ(pz*(1.0+dPtDivPt) - pT*dtheta/sin2theta);
 110
 111   Double_t mass2 = part->Mass()*part->Mass();
 112   Double_t e = mass2 + smearedPx*smearedPx +
 113                        smearedPy*smearedPy +
 114                        smearedPz*smearedPz;
 115
 116   smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(e));
 117 }
 118
 119 /****************************************************************/
 120 /****************************************************************/
 121 /****************************************************************/
 122
 123
 124 ClassImp(AliHBTCorrectQInvCorrelFctn)
 125
 126 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title):
 127   AliHBTOnePairFctn1D(name,title),
 128   fMeasCorrelFctn(0x0),
 129   fMeasNumer(0x0),
 130   fMeasDenom(0x0),
 131   fSmearedNumer(0x0),
 132   fSmearedDenom(0x0),
 133   fR2(0.0),
 134   fLambda(0.0),
 135   fRConvergenceTreshold(0.3),
 136   fLambdaConvergenceTreshold(0.05)
 137 {
 138
 139 }
 140 /******************************************************************/
 141
 142 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title):
 143   AliHBTOnePairFctn1D(name,title,
 144                       measqinv->GetNbinsX(),
 145           measqinv->GetXaxis()->GetXmax(),
 146           measqinv->GetXaxis()->GetXmin()),
 147   fMeasCorrelFctn(measqinv),
 148   fMeasNumer(0x0),
 149   fMeasDenom(0x0),
 150   fSmearedNumer(0x0),
 151   fSmearedDenom(0x0),
 152   fR2(0.0),
 153   fLambda(0.0),
 154   fRConvergenceTreshold(0.3),
 155   fLambdaConvergenceTreshold(0.05)
 156 {
 157
 158 }
 159 /******************************************************************/
 160
 161 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title,
 162                     Int_t nbins, Float_t maxXval, Float_t minXval):
 163   AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval),
 164   AliHBTCorrectedCorrelFctn(),
 165   fMeasCorrelFctn(0x0),
 166   fMeasNumer(0x0),
 167   fMeasDenom(0x0),
 168   fSmearedNumer(0x0),
 169   fSmearedDenom(0x0),
 170   fR2(0.0),
 171   fLambda(0.0),
 172   fRConvergenceTreshold(0.3),
 173   fLambdaConvergenceTreshold(0.05)
 174 {
 175 //ctor
 176 }
 177 /******************************************************************/
 178 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
 179   AliHBTOnePairFctn1D(in),
 180   AliHBTCorrectedCorrelFctn(),
 181   fMeasCorrelFctn(0x0),
 182   fMeasNumer(0x0),
 183   fMeasDenom(0x0),
 184   fSmearedNumer(0x0),
 185   fSmearedDenom(0x0),
 186   fR2(0.0),
 187   fLambda(0.0),
 188   fRConvergenceTreshold(0),
 189   fLambdaConvergenceTreshold(0)
 190 {
 191 //cpy ;ctor
 192  in.Copy(*this);
 193 }
 194
 195 AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
 196 {
 197 //dtor
 198   delete fMeasCorrelFctn;
 199   delete fSmearedNumer;
 200   delete fSmearedDenom;
 201   delete fMeasNumer;
 202   delete fMeasDenom;
 203 }
 204 /******************************************************************/
 205
 206 void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
 207 {
 208   AliHBTFunction1D::BuildHistos(nbins,max,min);
 209
 210   TString numstr = fName + " Smeared Numerator";  //title and name of the numerator histogram
 211   TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram
 212
 213   fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
 214   fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
 215   fSmearedNumer->Sumw2();
 216   fSmearedDenom->Sumw2();
 217
 218   if (fMeasCorrelFctn == 0x0)
 219    {
 220      numstr = fName + " Measured Numerator";  //title and name of the numerator histogram
 221      denstr = fName + " Measured Denominator";//title and name of the denominator histogram
 222
 223      fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
 224      fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
 225      fMeasNumer->Sumw2();
 226      fMeasDenom->Sumw2();
 227    }
 228 }
 229 /******************************************************************/
 230
 231 void AliHBTCorrectQInvCorrelFctn::Init()
 232 {
 233 //Init
 234   AliHBTOnePairFctn1D::Init();
 235   Info("Init","");
 236   fSmearedNumer->Reset();
 237   fSmearedDenom->Reset();
 238   if (fMeasNumer) fMeasNumer->Reset();
 239   if (fMeasDenom) fMeasDenom->Reset();
 240   fFittedR = -1.0;
 241   fFittedLambda = 0.0;
 242 }
 243 /******************************************************************/
 244
 245 void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
 246 {
 247  //Processes particles that originates from the same event
 248   if (fMeasNumer == 0x0) return;
 249   pair = CheckPair(pair);
 250   if( pair == 0x0) return;
 251   fMeasNumer->Fill(pair->GetQInv());
 252 }
 253 /******************************************************************/
 254
 255 void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
 256 {
 257 //Process different events
 258   static AliAODParticle part1, part2;
 259   static AliHBTPair smearedpair(&part1,&part2);
 260
 261   pair = CheckPair(pair);
 262   if( pair == 0x0) return;
 263   Double_t cc = GetCoulombCorrection(pair);
 264   Double_t qinv = pair->GetQInv();
 265
 266   //measured histogram -> if we are interested
 267   //only if fMeasCorrelFctn is not specified by user
 268   if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
 269
 270   Smear(pair,smearedpair);
 271   Double_t modelqinv = GetModelValue(qinv);
 272   //Ideal histogram
 273   fNumerator->Fill(qinv,modelqinv*cc);
 274   fDenominator->Fill(qinv,cc);
 275
 276   //Smeared histogram
 277   Double_t smearedqinv = smearedpair.GetQInv();
 278   fSmearedNumer->Fill(smearedqinv,modelqinv);
 279   Double_t smearedcc = GetCoulombCorrection(&smearedpair);
 280   fSmearedDenom->Fill(smearedqinv,smearedcc);
 281
 282 }
 283 /******************************************************************/
 284
 285 void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
 286 {
 287  //Sets Initial Values
 288   fLambda = lambda;
 289   fR2 = r*r;
 290 }
 291 /******************************************************************/
 292 void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
 293 {
 294   //makes measured correlation function
 295   delete fMeasCorrelFctn;
 296   fMeasCorrelFctn = 0x0;
 297
 298   if (fMeasNumer&&fMeasDenom)
 299    {
 300      Double_t measscale = Scale(fMeasNumer,fMeasDenom);
 301      if (measscale == 0.0)
 302       {
 303         Error("GetResult","GetRatio for measured CF returned 0.0");
 304         return;
 305       }
 306      TString str = fName + "measured ratio";
 307      fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
 308      fMeasCorrelFctn->SetTitle(str.Data());
 309      fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
 310    }
 311 }
 312
 313 TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
 314 {
 315   //In case we don't have yet Measured Correlation Function
 316   //Try to get it
 317   //result is
 318   //        N[meas]   N[ideal]/D[ideal]
 319   //C(Q) =  ------- * -----------------
 320   //        D[meas]   N[smear]/D[smear]
 321
 322   TString str;
 323   if (fMeasCorrelFctn == 0x0) MakeMeasCF();
 324
 325   if (fMeasCorrelFctn == 0x0)
 326    {
 327      Error("GetResult",
 328            "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
 329      return 0x0;
 330    }
 331
 332   TH1D* ideal = (TH1D*)GetRatio(Scale());
 333   if (ideal == 0x0)
 334    {
 335      Error("GetResult","Ratio of ideal histograms is null");
 336      return 0x0;
 337    }
 338   str = fName + " smeared ratio";
 339   TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
 340   smearedCF->SetTitle(str.Data());
 341   Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
 342   smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
 343
 344   str = fName + " product meas ideal CF";
 345   TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
 346   measideal->Multiply(ideal,fMeasCorrelFctn);
 347
 348   str = fName + " Corrected Result";
 349   TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
 350   result->SetTitle(str.Data());
 351
 352   Double_t resultscale = Scale(measideal,smearedCF);
 353   result->Divide(measideal,smearedCF,resultscale);
 354   return result;
 355
 356 }
 357 /******************************************************************/
 358
 359 void AliHBTCorrectQInvCorrelFctn::Fit()
 360 {
 361 //fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))
 362 //where [0] is lambda
 363 //      [1] is radius
 364 //   0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
 365
 366   Info("Fit","Before fFittedLambda = %f",fFittedLambda);
 367   Info("Fit","Before fFittedR = %f",fFittedR);
 368   TH1D* result = (TH1D*)GetResult();
 369   if (result == 0x0)
 370    {
 371      Error("Fit","Can not get result");
 372      return;
 373    }
 374   TF1* fitfctn = new TF1("fitfctn","1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))");
 375   fitfctn->SetParameter(0,1.0);
 376   fitfctn->SetParameter(1,6.0);
 377   Float_t max = result->GetXaxis()->GetXmax();
 378   Info("Fit","Max is %f",max);
 379   result->Fit(fitfctn,"0","",0.008,max);
 380   fFittedLambda = fitfctn->GetParameter(0);
 381   fFittedR = fitfctn->GetParameter(1);
 382   Info("Fit","After fFittedLambda = %f",fFittedLambda);
 383   Info("Fit","After fFittedR = %f",fFittedR);
 384   delete fitfctn;
 385   delete result;
 386 }
 387 /******************************************************************/
 388
 389 Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
 390 {
 391   //check if fitting was performed
 392   if (fFittedR <= 0.0)
 393    {
 394      Fit();//if not do fit first
 395      if (fFittedR <= 0.0)
 396       {
 397         Error("IsConverged","Fitting failed");
 398         return kFALSE;
 399       }
 400    }
 401
 402   Double_t guessedR = TMath::Sqrt(fR2);
 403
 404   Info("IsConverged","Fitted   lambda            : %8f Fitted  Radius             : %8f",fFittedLambda,fFittedR);
 405   Info("IsConverged","Guessed  lambda            : %8f Guessed Radius             : %8f",fLambda,guessedR);
 406   Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f",
 407        fLambdaConvergenceTreshold,fRConvergenceTreshold);
 408
 409   if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) &&
 410        (TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
 411     {
 412       Info("IsConverged","Cnvergence reached");
 413       return kTRUE;
 414     }
 415   else
 416    {
 417       Info("IsConverged","Cnvergence NOT reached");
 418       return kFALSE;
 419    }
 420 }
 421 /******************************************************************/
 422
 423 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
 424 {
 425  //Returns Fitted radius
 426   if (fFittedR <= 0.0) Fit();
 427   return fFittedR;
 428 }
 429 /******************************************************************/
 430
 431 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
 432 {
 433  //Returns Fitted Intercept paramter
 434   if (fFittedR <= 0.0) Fit();
 435   return fFittedLambda;
 436 }
 437 /******************************************************************/
 438
 439 void AliHBTCorrectQInvCorrelFctn::WriteAll()
 440 {
 441 //Writes function and all additional information
 442   Write();
 443   if (fMeasCorrelFctn) fMeasCorrelFctn->Write();
 444   if (fMeasNumer ) fMeasNumer->Write();
 445   if (fMeasDenom ) fMeasDenom->Write();
 446   if (fSmearedNumer) fSmearedNumer->Write();
 447   if (fSmearedDenom) fSmearedDenom->Write();
 448   if (fSmearedNumer && fSmearedDenom)
 449    {
 450     TString str = fName + " smeared ratio";
 451     TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
 452     smearedCF->SetTitle(str.Data());
 453     Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
 454     smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
 455    }
 456 }
 457