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://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 #include <AliAODParticle.h>
  41
  42
  43 ClassImp(AliHBTCorrectQInvCorrelFctn)
  44
  45 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title):
  46   AliHBTOnePairFctn1D(name,title),
  47   fMeasCorrelFctn(0x0),
  48   fMeasNumer(0x0),
  49   fMeasDenom(0x0),
  50   fSmearedNumer(0x0),
  51   fSmearedDenom(0x0),
  52   fDPtOverPtRMS(0.004),
  53   fThetaA(2.72e-03),
  54   fThetaB(1.87e-04),
  55   fThetaAlpha(-2.4),
  56   fPhiA(1.94e-03),
  57   fPhiB(2.5e-04),
  58   fPhiAlpha(-2.2),
  59   fR2(0.0),
  60   fLambda(0.0),
  61   fRConvergenceTreshold(0.3),
  62   fLambdaConvergenceTreshold(0.05)
  63 {
  64
  65 }
  66 /******************************************************************/
  67
  68 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title):
  69   AliHBTOnePairFctn1D(name,title,
  70                       measqinv->GetNbinsX(),
  71           measqinv->GetXaxis()->GetXmax(),
  72           measqinv->GetXaxis()->GetXmin()),
  73   fMeasCorrelFctn(measqinv),
  74   fMeasNumer(0x0),
  75   fMeasDenom(0x0),
  76   fSmearedNumer(0x0),
  77   fSmearedDenom(0x0),
  78   fDPtOverPtRMS(0.004),
  79   fThetaA(2.72e-03),
  80   fThetaB(1.87e-04),
  81   fThetaAlpha(-2.4),
  82   fPhiA(1.94e-03),
  83   fPhiB(2.5e-04),
  84   fPhiAlpha(-2.2),
  85   fR2(0.0),
  86   fLambda(0.0),
  87   fRConvergenceTreshold(0.3),
  88   fLambdaConvergenceTreshold(0.05)
  89 {
  90
  91 }
  92 /******************************************************************/
  93
  94 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title,
  95                     Int_t nbins, Float_t maxXval, Float_t minXval):
  96   AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval),
  97   fMeasCorrelFctn(0x0),
  98   fMeasNumer(0x0),
  99   fMeasDenom(0x0),
 100   fSmearedNumer(0x0),
 101   fSmearedDenom(0x0),
 102   fDPtOverPtRMS(0.004),
 103   fThetaA(2.72e-03),
 104   fThetaB(1.87e-04),
 105   fThetaAlpha(-2.4),
 106   fPhiA(1.94e-03),
 107   fPhiB(2.5e-04),
 108   fPhiAlpha(-2.2),
 109   fR2(0.0),
 110   fLambda(0.0),
 111   fRConvergenceTreshold(0.3),
 112   fLambdaConvergenceTreshold(0.05)
 113 {
 114 //ctor
 115 }
 116 /******************************************************************/
 117 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
 118   AliHBTOnePairFctn1D(in),
 119   fMeasCorrelFctn(0x0),
 120   fMeasNumer(0x0),
 121   fMeasDenom(0x0),
 122   fSmearedNumer(0x0),
 123   fSmearedDenom(0x0),
 124   fDPtOverPtRMS(0),
 125   fThetaA(0),
 126   fThetaB(0),
 127   fThetaAlpha(0),
 128   fPhiA(0),
 129   fPhiB(0),
 130   fPhiAlpha(0),
 131   fR2(0.0),
 132   fLambda(0.0),
 133   fRConvergenceTreshold(0),
 134   fLambdaConvergenceTreshold(0)
 135 {
 136 //cpy ;ctor
 137  in.Copy(*this);
 138 }
 139
 140 AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
 141 {
 142 //dtor
 143   delete fMeasCorrelFctn;
 144   delete fSmearedNumer;
 145   delete fSmearedDenom;
 146   delete fMeasNumer;
 147   delete fMeasDenom;
 148 }
 149 /******************************************************************/
 150
 151 void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
 152 {
 153   AliHBTFunction1D::BuildHistos(nbins,max,min);
 154
 155   TString numstr = fName + " Smeared Numerator";  //title and name of the numerator histogram
 156   TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram
 157
 158   fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
 159   fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
 160   fSmearedNumer->Sumw2();
 161   fSmearedDenom->Sumw2();
 162
 163   if (fMeasCorrelFctn == 0x0)
 164    {
 165      numstr = fName + " Measured Numerator";  //title and name of the numerator histogram
 166      denstr = fName + " Measured Denominator";//title and name of the denominator histogram
 167
 168      fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
 169      fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
 170      fMeasNumer->Sumw2();
 171      fMeasDenom->Sumw2();
 172    }
 173 }
 174 /******************************************************************/
 175
 176 void AliHBTCorrectQInvCorrelFctn::Init()
 177 {
 178 //Init
 179   AliHBTOnePairFctn1D::Init();
 180   Info("Init","");
 181   fSmearedNumer->Reset();
 182   fSmearedDenom->Reset();
 183   if (fMeasNumer) fMeasNumer->Reset();
 184   if (fMeasDenom) fMeasDenom->Reset();
 185   fFittedR = -1.0;
 186   fFittedLambda = 0.0;
 187 }
 188 /******************************************************************/
 189
 190 void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
 191 {
 192  //Processes particles that originates from the same event
 193   if (fMeasNumer == 0x0) return;
 194   pair = CheckPair(pair);
 195   if( pair == 0x0) return;
 196   fMeasNumer->Fill(pair->GetQInv());
 197 }
 198 /******************************************************************/
 199
 200 void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
 201 {
 202 //Process different events
 203   static AliAODParticle part1, part2;
 204   static AliHBTPair smearedpair(&part1,&part2);
 205
 206   pair = CheckPair(pair);
 207   if( pair == 0x0) return;
 208   Double_t cc = GetCoulombCorrection(pair);
 209   Double_t qinv = pair->GetQInv();
 210
 211   //measured histogram -> if we are interested
 212   //only if fMeasCorrelFctn is not specified by user
 213   if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
 214
 215   Smear(pair,smearedpair);
 216   Double_t modelqinv = GetModelValue(qinv);
 217   //Ideal histogram
 218   fNumerator->Fill(qinv,modelqinv*cc);
 219   fDenominator->Fill(qinv,cc);
 220
 221   //Smeared histogram
 222   Double_t smearedqinv = smearedpair.GetQInv();
 223   fSmearedNumer->Fill(smearedqinv,modelqinv);
 224   Double_t smearedcc = GetCoulombCorrection(&smearedpair);
 225   fSmearedDenom->Fill(smearedqinv,smearedcc);
 226
 227 }
 228 /******************************************************************/
 229 void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
 230 {
 231 //Smears pair
 232   Smear(pair->Particle1(),smeared.Particle1());
 233   Smear(pair->Particle2(),smeared.Particle2());
 234   smeared.Changed();
 235 }
 236 /******************************************************************/
 237
 238 void AliHBTCorrectQInvCorrelFctn::Smear(AliVAODParticle* part, AliVAODParticle* smeared)
 239 {
 240  //Smears momenta
 241   Double_t sin2theta = TMath::Sin(part->Theta());
 242   sin2theta = sin2theta*sin2theta;
 243   Double_t pt = part->Pt();
 244
 245   double dPtDivPt = gRandom->Gaus(0.0,fDPtOverPtRMS);
 246   double dphi = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fPhiAlpha));
 247   double dtheta = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fThetaAlpha));
 248
 249   Double_t smearedPx = part->Px()*(1.0+dPtDivPt) - part->Py()*dphi;
 250 //  fourmom.setX(px*(1.0+dPtDivPt) - py*dphi);
 251   Double_t smearedPy = part->Py()*(1.0+dPtDivPt) - part->Px()*dphi;
 252 //  fourmom.setY(py*(1.0+dPtDivPt) + px*dphi);
 253   Double_t smearedPz = part->Pz()*(1.0+dPtDivPt) - pt*dtheta/sin2theta;
 254 //  fourmom.setZ(pz*(1.0+dPtDivPt) - pT*dtheta/sin2theta);
 255
 256   Double_t mass2 = part->Mass()*part->Mass();
 257   Double_t e = mass2 + smearedPx*smearedPx +
 258                        smearedPy*smearedPy +
 259                        smearedPz*smearedPz;
 260
 261   smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(e));
 262 }
 263 /******************************************************************/
 264
 265 void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
 266 {
 267  //Sets Initial Values
 268   fLambda = lambda;
 269   fR2 = r*r;
 270 }
 271 /******************************************************************/
 272 void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
 273 {
 274   //makes measured correlation function
 275   delete fMeasCorrelFctn;
 276   fMeasCorrelFctn = 0x0;
 277
 278   if (fMeasNumer&&fMeasDenom)
 279    {
 280      Double_t measscale = Scale(fMeasNumer,fMeasDenom);
 281      if (measscale == 0.0)
 282       {
 283         Error("GetResult","GetRatio for measured CF returned 0.0");
 284         return;
 285       }
 286      TString str = fName + "measured ratio";
 287      fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
 288      fMeasCorrelFctn->SetTitle(str.Data());
 289      fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
 290    }
 291 }
 292
 293 TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
 294 {
 295   //In case we don't have yet Measured Correlation Function
 296   //Try to get it
 297   //result is
 298   //        N[meas]   N[ideal]/D[ideal]
 299   //C(Q) =  ------- * -----------------
 300   //        D[meas]   N[smear]/D[smear]
 301
 302   TString str;
 303   if (fMeasCorrelFctn == 0x0) MakeMeasCF();
 304
 305   if (fMeasCorrelFctn == 0x0)
 306    {
 307      Error("GetResult",
 308            "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
 309      return 0x0;
 310    }
 311
 312   TH1D* ideal = (TH1D*)GetRatio(Scale());
 313   if (ideal == 0x0)
 314    {
 315      Error("GetResult","Ratio of ideal histograms is null");
 316      return 0x0;
 317    }
 318   str = fName + " smeared ratio";
 319   TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
 320   smearedCF->SetTitle(str.Data());
 321   Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
 322   smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
 323
 324   str = fName + " product meas ideal CF";
 325   TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
 326   measideal->Multiply(ideal,fMeasCorrelFctn);
 327
 328   str = fName + " Corrected Result";
 329   TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
 330   result->SetTitle(str.Data());
 331
 332   Double_t resultscale = Scale(measideal,smearedCF);
 333   result->Divide(measideal,smearedCF,resultscale);
 334   return result;
 335
 336 }
 337 /******************************************************************/
 338
 339 void AliHBTCorrectQInvCorrelFctn::Fit()
 340 {
 341 //fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))
 342 //where [0] is lambda
 343 //      [1] is radius
 344 //   0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
 345
 346   Info("Fit","Before fFittedLambda = %f",fFittedLambda);
 347   Info("Fit","Before fFittedR = %f",fFittedR);
 348   TH1D* result = (TH1D*)GetResult();
 349   if (result == 0x0)
 350    {
 351      Error("Fit","Can not get result");
 352      return;
 353    }
 354   TF1* fitfctn = new TF1("fitfctn","1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))");
 355   fitfctn->SetParameter(0,1.0);
 356   fitfctn->SetParameter(1,6.0);
 357   Float_t max = result->GetXaxis()->GetXmax();
 358   Info("Fit","Max is %f",max);
 359   result->Fit(fitfctn,"0","",0.008,max);
 360   fFittedLambda = fitfctn->GetParameter(0);
 361   fFittedR = fitfctn->GetParameter(1);
 362   Info("Fit","After fFittedLambda = %f",fFittedLambda);
 363   Info("Fit","After fFittedR = %f",fFittedR);
 364   delete fitfctn;
 365   delete result;
 366 }
 367 /******************************************************************/
 368
 369 Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
 370 {
 371   //check if fitting was performed
 372   if (fFittedR <= 0.0)
 373    {
 374      Fit();//if not do fit first
 375      if (fFittedR <= 0.0)
 376       {
 377         Error("IsConverged","Fitting failed");
 378         return kFALSE;
 379       }
 380    }
 381
 382   Double_t guessedR = TMath::Sqrt(fR2);
 383
 384   Info("IsConverged","Fitted   lambda            : %8f Fitted  Radius             : %8f",fFittedLambda,fFittedR);
 385   Info("IsConverged","Guessed  lambda            : %8f Guessed Radius             : %8f",fLambda,guessedR);
 386   Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f",
 387        fLambdaConvergenceTreshold,fRConvergenceTreshold);
 388
 389   if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) &&
 390        (TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
 391     {
 392       Info("IsConverged","Cnvergence reached");
 393       return kTRUE;
 394     }
 395   else
 396    {
 397       Info("IsConverged","Cnvergence NOT reached");
 398       return kFALSE;
 399    }
 400 }
 401 /******************************************************************/
 402
 403 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
 404 {
 405  //Returns Fitted radius
 406   if (fFittedR <= 0.0) Fit();
 407   return fFittedR;
 408 }
 409 /******************************************************************/
 410
 411 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
 412 {
 413  //Returns Fitted Intercept paramter
 414   if (fFittedR <= 0.0) Fit();
 415   return fFittedLambda;
 416 }
 417 /******************************************************************/
 418
 419 void AliHBTCorrectQInvCorrelFctn::WriteAll()
 420 {
 421 //Writes function and all additional information
 422   Write();
 423   if (fMeasCorrelFctn) fMeasCorrelFctn->Write();
 424   if (fMeasNumer ) fMeasNumer->Write();
 425   if (fMeasDenom ) fMeasDenom->Write();
 426   if (fSmearedNumer) fSmearedNumer->Write();
 427   if (fSmearedDenom) fSmearedDenom->Write();
 428   if (fSmearedNumer && fSmearedDenom)
 429    {
 430     TString str = fName + " smeared ratio";
 431     TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
 432     smearedCF->SetTitle(str.Data());
 433     Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
 434     smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
 435    }
 436 }
 437