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   fFittedR(0),
 136   fFittedLambda(0),
 137   fRConvergenceTreshold(0.3),
 138   fLambdaConvergenceTreshold(0.05)
 139 {
 140
 141 }
 142 /******************************************************************/
 143
 144 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title):
 145   AliHBTOnePairFctn1D(name,title,
 146                       measqinv->GetNbinsX(),
 147           measqinv->GetXaxis()->GetXmax(),
 148           measqinv->GetXaxis()->GetXmin()),
 149   fMeasCorrelFctn(measqinv),
 150   fMeasNumer(0x0),
 151   fMeasDenom(0x0),
 152   fSmearedNumer(0x0),
 153   fSmearedDenom(0x0),
 154   fR2(0.0),
 155   fLambda(0.0),
 156   fFittedR(0),
 157   fFittedLambda(0),
 158   fRConvergenceTreshold(0.3),
 159   fLambdaConvergenceTreshold(0.05)
 160 {
 161
 162 }
 163 /******************************************************************/
 164
 165 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title,
 166                     Int_t nbins, Float_t maxXval, Float_t minXval):
 167   AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval),
 168   AliHBTCorrectedCorrelFctn(),
 169   fMeasCorrelFctn(0x0),
 170   fMeasNumer(0x0),
 171   fMeasDenom(0x0),
 172   fSmearedNumer(0x0),
 173   fSmearedDenom(0x0),
 174   fR2(0.0),
 175   fLambda(0.0),
 176   fFittedR(0),
 177   fFittedLambda(0),
 178   fRConvergenceTreshold(0.3),
 179   fLambdaConvergenceTreshold(0.05)
 180 {
 181 //ctor
 182 }
 183 /******************************************************************/
 184 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
 185   AliHBTOnePairFctn1D(in),
 186   AliHBTCorrectedCorrelFctn(),
 187   fMeasCorrelFctn(0x0),
 188   fMeasNumer(0x0),
 189   fMeasDenom(0x0),
 190   fSmearedNumer(0x0),
 191   fSmearedDenom(0x0),
 192   fR2(0.0),
 193   fLambda(0.0),
 194   fFittedR(0),
 195   fFittedLambda(0),
 196   fRConvergenceTreshold(0),
 197   fLambdaConvergenceTreshold(0)
 198 {
 199 //cpy ;ctor
 200  in.Copy(*this);
 201 }
 202
 203 AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
 204 {
 205 //dtor
 206   delete fMeasCorrelFctn;
 207   delete fSmearedNumer;
 208   delete fSmearedDenom;
 209   delete fMeasNumer;
 210   delete fMeasDenom;
 211 }
 212 /******************************************************************/
 213
 214 void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
 215 {
 216   AliHBTFunction1D::BuildHistos(nbins,max,min);
 217
 218   TString numstr = fName + " Smeared Numerator";  //title and name of the numerator histogram
 219   TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram
 220
 221   fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
 222   fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
 223   fSmearedNumer->Sumw2();
 224   fSmearedDenom->Sumw2();
 225
 226   if (fMeasCorrelFctn == 0x0)
 227    {
 228      numstr = fName + " Measured Numerator";  //title and name of the numerator histogram
 229      denstr = fName + " Measured Denominator";//title and name of the denominator histogram
 230
 231      fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
 232      fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
 233      fMeasNumer->Sumw2();
 234      fMeasDenom->Sumw2();
 235    }
 236 }
 237 /******************************************************************/
 238
 239 void AliHBTCorrectQInvCorrelFctn::Init()
 240 {
 241 //Init
 242   AliHBTOnePairFctn1D::Init();
 243   Info("Init","");
 244   fSmearedNumer->Reset();
 245   fSmearedDenom->Reset();
 246   if (fMeasNumer) fMeasNumer->Reset();
 247   if (fMeasDenom) fMeasDenom->Reset();
 248   fFittedR = -1.0;
 249   fFittedLambda = 0.0;
 250 }
 251 /******************************************************************/
 252
 253 void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
 254 {
 255  //Processes particles that originates from the same event
 256   if (fMeasNumer == 0x0) return;
 257   pair = CheckPair(pair);
 258   if( pair == 0x0) return;
 259   fMeasNumer->Fill(pair->GetQInv());
 260 }
 261 /******************************************************************/
 262
 263 void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
 264 {
 265 //Process different events
 266   static AliAODParticle part1, part2;
 267   static AliHBTPair smearedpair(&part1,&part2);
 268
 269   pair = CheckPair(pair);
 270   if( pair == 0x0) return;
 271   Double_t cc = GetCoulombCorrection(pair);
 272   Double_t qinv = pair->GetQInv();
 273
 274   //measured histogram -> if we are interested
 275   //only if fMeasCorrelFctn is not specified by user
 276   if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
 277
 278   Smear(pair,smearedpair);
 279   Double_t modelqinv = GetModelValue(qinv);
 280   //Ideal histogram
 281   fNumerator->Fill(qinv,modelqinv*cc);
 282   fDenominator->Fill(qinv,cc);
 283
 284   //Smeared histogram
 285   Double_t smearedqinv = smearedpair.GetQInv();
 286   fSmearedNumer->Fill(smearedqinv,modelqinv);
 287   Double_t smearedcc = GetCoulombCorrection(&smearedpair);
 288   fSmearedDenom->Fill(smearedqinv,smearedcc);
 289
 290 }
 291 /******************************************************************/
 292
 293 void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
 294 {
 295  //Sets Initial Values
 296   fLambda = lambda;
 297   fR2 = r*r;
 298 }
 299 /******************************************************************/
 300 void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
 301 {
 302   //makes measured correlation function
 303   delete fMeasCorrelFctn;
 304   fMeasCorrelFctn = 0x0;
 305
 306   if (fMeasNumer&&fMeasDenom)
 307    {
 308      Double_t measscale = Scale(fMeasNumer,fMeasDenom);
 309      if (measscale == 0.0)
 310       {
 311         Error("GetResult","GetRatio for measured CF returned 0.0");
 312         return;
 313       }
 314      TString str = fName + "measured ratio";
 315      fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
 316      fMeasCorrelFctn->SetTitle(str.Data());
 317      fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
 318    }
 319 }
 320
 321 TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
 322 {
 323   //In case we don't have yet Measured Correlation Function
 324   //Try to get it
 325   //result is
 326   //        N[meas]   N[ideal]/D[ideal]
 327   //C(Q) =  ------- * -----------------
 328   //        D[meas]   N[smear]/D[smear]
 329
 330   TString str;
 331   if (fMeasCorrelFctn == 0x0) MakeMeasCF();
 332
 333   if (fMeasCorrelFctn == 0x0)
 334    {
 335      Error("GetResult",
 336            "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
 337      return 0x0;
 338    }
 339
 340   TH1D* ideal = (TH1D*)GetRatio(Scale());
 341   if (ideal == 0x0)
 342    {
 343      Error("GetResult","Ratio of ideal histograms is null");
 344      return 0x0;
 345    }
 346   str = fName + " smeared ratio";
 347   TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
 348   smearedCF->SetTitle(str.Data());
 349   Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
 350   smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
 351
 352   str = fName + " product meas ideal CF";
 353   TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
 354   measideal->Multiply(ideal,fMeasCorrelFctn);
 355
 356   str = fName + " Corrected Result";
 357   TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
 358   result->SetTitle(str.Data());
 359
 360   Double_t resultscale = Scale(measideal,smearedCF);
 361   result->Divide(measideal,smearedCF,resultscale);
 362   return result;
 363
 364 }
 365 /******************************************************************/
 366
 367 void AliHBTCorrectQInvCorrelFctn::Fit()
 368 {
 369 //fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))
 370 //where [0] is lambda
 371 //      [1] is radius
 372 //   0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
 373
 374   Info("Fit","Before fFittedLambda = %f",fFittedLambda);
 375   Info("Fit","Before fFittedR = %f",fFittedR);
 376   TH1D* result = (TH1D*)GetResult();
 377   if (result == 0x0)
 378    {
 379      Error("Fit","Can not get result");
 380      return;
 381    }
 382   TF1* fitfctn = new TF1("fitfctn","1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))");
 383   fitfctn->SetParameter(0,1.0);
 384   fitfctn->SetParameter(1,6.0);
 385   Float_t max = result->GetXaxis()->GetXmax();
 386   Info("Fit","Max is %f",max);
 387   result->Fit(fitfctn,"0","",0.008,max);
 388   fFittedLambda = fitfctn->GetParameter(0);
 389   fFittedR = fitfctn->GetParameter(1);
 390   Info("Fit","After fFittedLambda = %f",fFittedLambda);
 391   Info("Fit","After fFittedR = %f",fFittedR);
 392   delete fitfctn;
 393   delete result;
 394 }
 395 /******************************************************************/
 396
 397 Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
 398 {
 399   //check if fitting was performed
 400   if (fFittedR <= 0.0)
 401    {
 402      Fit();//if not do fit first
 403      if (fFittedR <= 0.0)
 404       {
 405         Error("IsConverged","Fitting failed");
 406         return kFALSE;
 407       }
 408    }
 409
 410   Double_t guessedR = TMath::Sqrt(fR2);
 411
 412   Info("IsConverged","Fitted   lambda            : %8f Fitted  Radius             : %8f",fFittedLambda,fFittedR);
 413   Info("IsConverged","Guessed  lambda            : %8f Guessed Radius             : %8f",fLambda,guessedR);
 414   Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f",
 415        fLambdaConvergenceTreshold,fRConvergenceTreshold);
 416
 417   if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) &&
 418        (TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
 419     {
 420       Info("IsConverged","Cnvergence reached");
 421       return kTRUE;
 422     }
 423   else
 424    {
 425       Info("IsConverged","Cnvergence NOT reached");
 426       return kFALSE;
 427    }
 428 }
 429 /******************************************************************/
 430
 431 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
 432 {
 433  //Returns Fitted radius
 434   if (fFittedR <= 0.0) Fit();
 435   return fFittedR;
 436 }
 437 /******************************************************************/
 438
 439 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
 440 {
 441  //Returns Fitted Intercept paramter
 442   if (fFittedR <= 0.0) Fit();
 443   return fFittedLambda;
 444 }
 445 /******************************************************************/
 446
 447 void AliHBTCorrectQInvCorrelFctn::WriteAll()
 448 {
 449 //Writes function and all additional information
 450   Write();
 451   if (fMeasCorrelFctn) fMeasCorrelFctn->Write();
 452   if (fMeasNumer ) fMeasNumer->Write();
 453   if (fMeasDenom ) fMeasDenom->Write();
 454   if (fSmearedNumer) fSmearedNumer->Write();
 455   if (fSmearedDenom) fSmearedDenom->Write();
 456   if (fSmearedNumer && fSmearedDenom)
 457    {
 458     TString str = fName + " smeared ratio";
 459     TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
 460     smearedCF->SetTitle(str.Data());
 461     Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
 462     smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
 463    }
 464 }
 465