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