From 83b33650b6b608a387ad8de7dcb6065e745fed9e Mon Sep 17 00:00:00 2001 From: skowron Date: Sun, 7 Sep 2003 17:20:10 +0000 Subject: [PATCH] Classes for Corrections with Lisa algorithm added --- HBTAN/AliHBTCorrectCorrelFctn.cxx | 441 ++++++++++++++++++++++++++++++ HBTAN/AliHBTCorrectCorrelFctn.h | 132 +++++++++ HBTAN/HBTAnalysisLinkDef.h | 3 + HBTAN/hbtcorrections.C | 81 ++++++ HBTAN/libHBTAN.pkg | 2 +- 5 files changed, 658 insertions(+), 1 deletion(-) create mode 100644 HBTAN/AliHBTCorrectCorrelFctn.cxx create mode 100644 HBTAN/AliHBTCorrectCorrelFctn.h create mode 100644 HBTAN/hbtcorrections.C diff --git a/HBTAN/AliHBTCorrectCorrelFctn.cxx b/HBTAN/AliHBTCorrectCorrelFctn.cxx new file mode 100644 index 00000000000..c50174c0a57 --- /dev/null +++ b/HBTAN/AliHBTCorrectCorrelFctn.cxx @@ -0,0 +1,441 @@ +#include "AliHBTCorrectCorrelFctn.h" +//____________________ +/////////////////////////////////////////////////////// +// // +// AliHBTCorrectQInvCorrelFctn // +// // +// Class for calculating Q Invariant correlation // +// taking to the account resolution of the // +// detector and coulomb effects. // +// Implemented on basis of STAR procedure // +// implemented and described by Michael A. Lisa // +// // +// Piotr.Skowronski@cern.ch // +// http://alisoft.cern.ch/people/skowron/analyzer // +// // +/////////////////////////////////////////////////////// + +//Parameters fit from pi+ pi+ resolution analysis for pair with qinv<50MeV +// chi2/NFD = 97/157 +// FCN=98.0971 FROM MIGRAD STATUS=CONVERGED 332 CALLS 333 TOTAL +// EDM=2.13364e-12 STRATEGY= 1 ERROR MATRIX UNCERTAINTY 2.4 per cent +// EXT PARAMETER STEP FIRST +// NO. NAME VALUE ERROR SIZE DERIVATIVE +// 1 fThetaA 2.72546e-03 1.43905e-04 -0.00000e+00 5.54375e-02 +// 2 fThetaB 1.87116e-04 5.11862e-05 0.00000e+00 3.66500e-01 +// 3 fThetaAlpha -2.36868e+00 1.83230e-01 0.00000e+00 -7.01301e-05 + +// FCN=120.603 FROM MIGRAD STATUS=CONVERGED 117 CALLS 118 TOTAL +// EDM=2.5153e-12 STRATEGY= 1 ERROR MATRIX UNCERTAINTY 2.4 per cent +// EXT PARAMETER STEP FIRST +// NO. NAME VALUE ERROR SIZE DERIVATIVE +// 1 fPhiA 1.93913e-03 1.31059e-04 -0.00000e+00 -5.87280e-02 +// 2 fPhiA 2.48687e-04 5.41251e-05 0.00000e+00 -1.87361e-01 +// 3 fPhiAlpha -2.22649e+00 1.44503e-01 0.00000e+00 8.97538e-06 + +#include +#include +#include +#include +ClassImp(AliHBTCorrectQInvCorrelFctn) + +AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title): + AliHBTOnePairFctn1D(name,title), + fMeasCorrelFctn(0x0), + fMeasNumer(0x0), + fMeasDenom(0x0), + fSmearedNumer(0x0), + fSmearedDenom(0x0), + fDPtOverPtRMS(0.004), + fThetaA(2.72e-03), + fThetaB(1.87e-04), + fThetaAlpha(-2.4), + fPhiA(1.94e-03), + fPhiB(2.5e-04), + fPhiAlpha(-2.2), + fR2(0.0), + fLambda(0.0), + fRConvergenceTreshold(0.3), + fLambdaConvergenceTreshold(0.05) +{ + +} +/******************************************************************/ + +AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title): + AliHBTOnePairFctn1D(name,title, + measqinv->GetNbinsX(), + measqinv->GetXaxis()->GetXmax(), + measqinv->GetXaxis()->GetXmin()), + fMeasCorrelFctn(measqinv), + fMeasNumer(0x0), + fMeasDenom(0x0), + fSmearedNumer(0x0), + fSmearedDenom(0x0), + fDPtOverPtRMS(0.004), + fThetaA(2.72e-03), + fThetaB(1.87e-04), + fThetaAlpha(-2.4), + fPhiA(1.94e-03), + fPhiB(2.5e-04), + fPhiAlpha(-2.2), + fR2(0.0), + fLambda(0.0), + fRConvergenceTreshold(0.3), + fLambdaConvergenceTreshold(0.05) +{ + +} +/******************************************************************/ + +AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title, + Int_t nbins, Float_t maxXval, Float_t minXval): + AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval), + fMeasCorrelFctn(0x0), + fMeasNumer(0x0), + fMeasDenom(0x0), + fSmearedNumer(0x0), + fSmearedDenom(0x0), + fDPtOverPtRMS(0.004), + fThetaA(2.72e-03), + fThetaB(1.87e-04), + fThetaAlpha(-2.4), + fPhiA(1.94e-03), + fPhiB(2.5e-04), + fPhiAlpha(-2.2), + fR2(0.0), + fLambda(0.0), + fRConvergenceTreshold(0.3), + fLambdaConvergenceTreshold(0.05) +{ +} +/******************************************************************/ + +AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn() +{ + delete fMeasCorrelFctn; + delete fSmearedNumer; + delete fSmearedDenom; + delete fMeasNumer; + delete fMeasDenom; +} +/******************************************************************/ + +void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min) +{ + AliHBTFunction1D::BuildHistos(nbins,max,min); + + TString numstr = fName + " Smeared Numerator"; //title and name of the numerator histogram + TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram + + fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max); + fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max); + fSmearedNumer->Sumw2(); + fSmearedDenom->Sumw2(); + + if (fMeasCorrelFctn == 0x0) + { + numstr = fName + " Measured Numerator"; //title and name of the numerator histogram + denstr = fName + " Measured Denominator";//title and name of the denominator histogram + + fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max); + fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max); + fMeasNumer->Sumw2(); + fMeasDenom->Sumw2(); + } +} +/******************************************************************/ + +void AliHBTCorrectQInvCorrelFctn::Init() +{ + AliHBTOnePairFctn1D::Init(); + Info("Init",""); + fSmearedNumer->Reset(); + fSmearedDenom->Reset(); + if (fMeasNumer) fMeasNumer->Reset(); + if (fMeasDenom) fMeasDenom->Reset(); + fFittedR = -1.0; + fFittedLambda = 0.0; +} +/******************************************************************/ + +void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair) +{ + if (fMeasNumer == 0x0) return; + pair = CheckPair(pair); + if( pair == 0x0) return; + fMeasNumer->Fill(pair->GetQInv()); +} +/******************************************************************/ + +void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair) +{ +//Process different events + static AliHBTParticle part1, part2; + static AliHBTPair smearedpair(&part1,&part2); + + pair = CheckPair(pair); + if( pair == 0x0) return; + Double_t cc = GetCoulombCorrection(pair); + Double_t qinv = pair->GetQInv(); + + //measured histogram -> if we are interested + //only if fMeasCorrelFctn is not specified by user + if (fMeasDenom) fMeasDenom->Fill(qinv,cc); + + Smear(pair,smearedpair); + Double_t modelqinv = GetModelValue(qinv); + //Ideal histogram + fNumerator->Fill(qinv,modelqinv*cc); + fDenominator->Fill(qinv,cc); + + //Smeared histogram + Double_t smearedqinv = smearedpair.GetQInv(); + fSmearedNumer->Fill(smearedqinv,modelqinv); + Double_t smearedcc = GetCoulombCorrection(&smearedpair); + fSmearedDenom->Fill(smearedqinv,smearedcc); + +} +/******************************************************************/ +void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared) +{ +//Smears pair + Smear(pair->Particle1(),smeared.Particle1()); + Smear(pair->Particle2(),smeared.Particle2()); + smeared.Changed(); +} +/******************************************************************/ + +void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTParticle* part, AliHBTParticle* smeared) +{ + Double_t sin2theta = TMath::Sin(part->Theta()); + sin2theta = sin2theta*sin2theta; + Double_t pt = part->Pt(); + + double DpT_div_pT = gRandom->Gaus(0.0,fDPtOverPtRMS); + double Dphi = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fPhiAlpha)); + double Dtheta = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fThetaAlpha)); + + Double_t smearedPx = part->Px()*(1.0+DpT_div_pT) - part->Py()*Dphi; +// fourmom.setX(px*(1.0+DpT_div_pT) - py*Dphi); + Double_t smearedPy = part->Py()*(1.0+DpT_div_pT) - part->Px()*Dphi; +// fourmom.setY(py*(1.0+DpT_div_pT) + px*Dphi); + Double_t smearedPz = part->Pz()*(1.0+DpT_div_pT) - pt*Dtheta/sin2theta; +// fourmom.setZ(pz*(1.0+DpT_div_pT) - pT*Dtheta/sin2theta); + + Double_t mass2 = part->GetMass()*part->GetMass(); + Double_t E = mass2 + smearedPx*smearedPx + + smearedPy*smearedPy + + smearedPz*smearedPz; + + smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(E)); +} +/******************************************************************/ + +void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r) +{ + fLambda = lambda; + fR2 = r*r; +} +/******************************************************************/ +void AliHBTCorrectQInvCorrelFctn::MakeMeasCF() +{ + //makes measured correlation function + delete fMeasCorrelFctn; + fMeasCorrelFctn = 0x0; + + if (fMeasNumer&&fMeasDenom) + { + Double_t measscale = Scale(fMeasNumer,fMeasDenom); + if (measscale == 0.0) + { + Error("GetResult","GetRatio for measured CF returned 0.0"); + return; + } + TString str = fName + "measured ratio"; + fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data()); + fMeasCorrelFctn->SetTitle(str.Data()); + fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale); + } +} + +TH1* AliHBTCorrectQInvCorrelFctn::GetResult() +{ + //In case we don't have yet Measured Correlation Function + //Try to get it + //result is + // N[meas] N[ideal]/D[ideal] + //C(Q) = ------- * ----------------- + // D[meas] N[smear]/D[smear] + + TString str; + if (fMeasCorrelFctn == 0x0) MakeMeasCF(); + + if (fMeasCorrelFctn == 0x0) + { + Error("GetResult", + "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null"); + return 0x0; + } + + TH1D* ideal = (TH1D*)GetRatio(Scale()); + if (ideal == 0x0) + { + Error("GetResult","Ratio of ideal histograms is null"); + return 0x0; + } + str = fName + " smeared ratio"; + TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data()); + smearedCF->SetTitle(str.Data()); + Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom); + smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale); + + str = fName + " product meas ideal CF"; + TH1D* measideal = (TH1D*)ideal->Clone(str.Data()); + measideal->Multiply(ideal,fMeasCorrelFctn); + + str = fName + " Corrected Result"; + TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data()); + result->SetTitle(str.Data()); + + Double_t resultscale = Scale(measideal,smearedCF); + result->Divide(measideal,smearedCF,resultscale); + return result; + +} +/******************************************************************/ + +void AliHBTCorrectQInvCorrelFctn::Fit() +{ +//fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329)) +//where [0] is lambda +// [1] is radius +// 0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI + + Info("Fit","Before fFittedLambda = %f",fFittedLambda); + Info("Fit","Before fFittedR = %f",fFittedR); + TH1D* result = (TH1D*)GetResult(); + if (result == 0x0) + { + Error("Fit","Can not get result"); + return; + } + TF1* fitfctn = new TF1("fitfctn","1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))"); + fitfctn->SetParameter(0,1.0); + fitfctn->SetParameter(1,6.0); + Float_t max = result->GetXaxis()->GetXmax(); + Info("Fit","Max is %f",max); + result->Fit(fitfctn,"0","",0.008,max); + fFittedLambda = fitfctn->GetParameter(0); + fFittedR = fitfctn->GetParameter(1); + Info("Fit","After fFittedLambda = %f",fFittedLambda); + Info("Fit","After fFittedR = %f",fFittedR); + delete fitfctn; + delete result; +} +/******************************************************************/ + +Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged() +{ + //check if fitting was performed + if (fFittedR <= 0.0) + { + Fit();//if not do fit first + if (fFittedR <= 0.0) + { + Error("IsConverged","Fitting failed"); + return kFALSE; + } + } + + Double_t guessedR = TMath::Sqrt(fR2); + + Info("IsConverged","Fitted lambda : %8f Fitted Radius : %8f",fFittedLambda,fFittedR); + Info("IsConverged","Guessed lambda : %8f Guessed Radius : %8f",fLambda,guessedR); + Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f", + fLambdaConvergenceTreshold,fRConvergenceTreshold); + + if ( (TMath::Abs(fLambda-fFittedLambda)Write(); + if (fMeasNumer ) fMeasNumer->Write(); + if (fMeasDenom ) fMeasDenom->Write(); + if (fSmearedNumer) fSmearedNumer->Write(); + if (fSmearedDenom) fSmearedDenom->Write(); + if (fSmearedNumer && fSmearedDenom) + { + TString str = fName + " smeared ratio"; + TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data()); + smearedCF->SetTitle(str.Data()); + Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom); + smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale); + } +} + +/******************************************************************/ +/******************************************************************/ +/******************************************************************/ + +//____________________ +/////////////////////////////////////////////////////// +// // +// AliHBTCorrectQ3DCorrelFctn // +// // +// Class for calculating Q Invariant correlation // +// taking to the account resolution of the // +// detector and coulomb effects. // +// // +/////////////////////////////////////////////////////// + + +AliHBTCorrectQ3DCorrelFctn::AliHBTCorrectQ3DCorrelFctn(const char* name, const char* title): + AliHBTOnePairFctn3D(name,title), + fMeasCorrelFctn(0x0), + fSmearedNumer(0x0), + fSmearedDenom(0x0), + fMeasNumer(0x0), + fMeasDenom(0x0) +{ + +} +/******************************************************************/ + +AliHBTCorrectQ3DCorrelFctn::~AliHBTCorrectQ3DCorrelFctn() +{ + delete fMeasCorrelFctn; + delete fSmearedNumer; + delete fSmearedDenom; + delete fMeasNumer; + delete fMeasDenom; +} + +/******************************************************************/ diff --git a/HBTAN/AliHBTCorrectCorrelFctn.h b/HBTAN/AliHBTCorrectCorrelFctn.h new file mode 100644 index 00000000000..01f0eb79730 --- /dev/null +++ b/HBTAN/AliHBTCorrectCorrelFctn.h @@ -0,0 +1,132 @@ +#ifndef ALIHBTCORRECTCORRELFCTN_H +#define ALIHBTCORRECTCORRELFCTN_H +//____________________ +/////////////////////////////////////////////////////// +// // +// AliHBTCorrectQInvCorrelFctn // +// // +// Class for calculating Q Invariant correlation // +// taking to the account resolution of the // +// detector and coulomb effects. // +// // +/////////////////////////////////////////////////////// + +#include "AliHBTFunction.h" + + class AliHBTCorrectQInvCorrelFctn: public AliHBTOnePairFctn1D +{ + public: + AliHBTCorrectQInvCorrelFctn(const char* name = "qinvcorrectedCF", + const char* title= "Corrected Q_{inv} Correlation Fonction"); + + AliHBTCorrectQInvCorrelFctn(const char* name, const char* title, + Int_t nbins, Float_t maxXval, Float_t minXval); + + AliHBTCorrectQInvCorrelFctn(TH1D* measqinv, + const char* name = "qinvcorrectedCF", + const char* title= "Corrected Q_{inv} Correlation Fonction"); + virtual ~AliHBTCorrectQInvCorrelFctn(); + void SetInitialValues(Double_t lambda, Double_t r); + void Init(); + void ProcessSameEventParticles(AliHBTPair* pair);//process particles from same event (real pair) + void ProcessDiffEventParticles(AliHBTPair* pair);//process particles coming from different events (mixed pairs) + void SetMeasuredHistogram(TH1D* meas){fMeasCorrelFctn = meas;} + TH1* GetResult();//returns the result histogram + Double_t GetRadius()const{ return TMath::Sqrt(fR2);}//returns assumed radius + Double_t GetLambda()const{ return fLambda;}//retutrns assumed intercept parameter + void SetRadiusConvergenceTreshold(Double_t ct){fRConvergenceTreshold=ct;}//if fitted and assumed R us different less then that number con + void SetLambdaConvergenceTreshold(Double_t ct){fLambdaConvergenceTreshold=ct;} + Bool_t IsConverged(); + void Fit(); + Double_t GetFittedRadius(); + Double_t GetFittedLambda(); + void WriteAll(); + void SetMeasNum(TH1D* measnum){fMeasNumer = measnum;} + void SetMeasDen(TH1D* h){fMeasDenom = h;} + void MakeMeasCF(); + protected: + virtual void BuildHistos(Int_t nbins, Float_t max, Float_t min); + Double_t GetCoulombCorrection(AliHBTPair* pair){return 1.0;} + Double_t GetValue(AliHBTPair * pair){return pair->GetQInv();} + void Smear(AliHBTPair* pair,AliHBTPair& smeared); + void Smear(AliHBTParticle* part, AliHBTParticle* smeared); + Double_t GetModelValue(Double_t qinv); + + //Our ideal numerator + TH1D* fMeasCorrelFctn; + TH1D* fMeasNumer; + TH1D* fMeasDenom; + + TH1D* fSmearedNumer; //! Numerator of smeard q + TH1D* fSmearedDenom; //! Denominator of smeard q + + //Parameters of Pt RMS + //linear dependence dPt/Pt from Pt itself + Float_t fDPtOverPtRMS; //RMS of dPt/Pt + + //We assume that RMS of Theta and Phisangle depends on Pt Like A+B*(Pt)^Alpha + //Idea copied from Star HBT Maker (Fabrice Retiere) + //Parameters comes from Monte Carlo Resolution Analysis + + Float_t fThetaA; //"A" parameter of theta RMS dependence + Float_t fThetaB; //"B" parameter of theta RMS dependence + Float_t fThetaAlpha; //"Alpha" parameter (power) of theta RMS dependence + + Float_t fPhiA;//"A" parameter of phi RMS dependence + Float_t fPhiB;//"B" parameter of phi RMS dependence + Float_t fPhiAlpha;//"Alpha" parameter (power) of phi RMS dependence + + Double_t fR2;//square of radius + Double_t fLambda;//Interception parameter + + Double_t fFittedR;//fitted radius + Double_t fFittedLambda;//fitted Interception parameter + + Float_t fRConvergenceTreshold; + Float_t fLambdaConvergenceTreshold; + private: + public: + ClassDef(AliHBTCorrectQInvCorrelFctn,1) +}; + +inline +Double_t AliHBTCorrectQInvCorrelFctn::GetModelValue(Double_t qinv) +{ + //factor 0.038936366329 conected with units change GeV<->SI + return 1.0 + fLambda*TMath::Exp(-fR2*qinv*qinv/0.038936366329); +} + +//____________________ +/////////////////////////////////////////////////////// +// // +// AliHBTCorrectQ3DCorrelFctn // +// // +// Class for calculating Q Invariant correlation // +// taking to the account resolution of the // +// detector and coulomb effects. // +// // +/////////////////////////////////////////////////////// + +class AliHBTCorrectQ3DCorrelFctn: public AliHBTOnePairFctn3D +{ + public: + AliHBTCorrectQ3DCorrelFctn(const char* name = "qinvcorrectedCF", + const char* title= "Corrected Q_{inv} Correlation Fonction"); + virtual ~AliHBTCorrectQ3DCorrelFctn(); + + protected: + TH3D* fMeasCorrelFctn; + + TH3D* fSmearedNumer; //! Numerator of smeard q + TH3D* fSmearedDenom; //! Denominator of smeard q + TH3D* fMeasNumer; //! Numerator of ideal q calculated on basis of model equation + TH3D* fMeasDenom; //! Denominator of ideal q calculated on basis of model equation + + + private: + + public: + ClassDef(AliHBTCorrectQ3DCorrelFctn,1) +}; + +#endif diff --git a/HBTAN/HBTAnalysisLinkDef.h b/HBTAN/HBTAnalysisLinkDef.h index 4c8c7c3fc95..c0efe59b2d8 100644 --- a/HBTAN/HBTAnalysisLinkDef.h +++ b/HBTAN/HBTAnalysisLinkDef.h @@ -85,6 +85,9 @@ #pragma link C++ class AliHBTInvMassCorrelFctn+; #pragma link C++ class AliHBTCorrFitFctn+; +#pragma link C++ class AliHBTCorrectQInvCorrelFctn+; +#pragma link C++ class AliHBTCorrectQ3DCorrelFctn+; + #pragma link C++ class AliHBTKtResolVsQInvFctn+; #pragma link C++ class AliHBTQInvResolVsQInvFctn+; diff --git a/HBTAN/hbtcorrections.C b/HBTAN/hbtcorrections.C new file mode 100644 index 00000000000..bd9ebf88a31 --- /dev/null +++ b/HBTAN/hbtcorrections.C @@ -0,0 +1,81 @@ +void hbtcorrections(Int_t first = -1,Int_t last = -1) +{ + + const char* basedir="."; + const char* serie=""; + const char* field = ""; + + AliHBTReader* reader = new AliHBTReaderInternal("ESD.old.root"); + TObjArray* dirs=0; + if ( (first >= 0) && (last>=0) && ( (last-first)>=0 ) ) + {//read from many dirs dirs + char buff[50]; + dirs = new TObjArray(last-first+1); + for (Int_t i = first; i<=last; i++) + { + sprintf(buff,"%s/%s/%s/%d",basedir,field,serie,i); + TObjString *odir= new TObjString(buff); + dirs->Add(odir); + } + } + reader->SetDirs(dirs); + + AliHBTParticleCut* readerpartcut= new AliHBTParticleCut(); + readerpartcut->SetPtRange(0.0,1.0); + readerpartcut->SetPID(kPiPlus); + reader->AddParticleCut(readerpartcut);//read this particle type with this cut + + AliHBTAnalysis* analysis = new AliHBTAnalysis(); + analysis->SetReader(reader); + analysis->SetDisplayInfo(100000); + AliHBTPairCut *paircut = new AliHBTPairCut(); + paircut->SetQInvRange(0.0,0.15); + analysis->SetGlobalPairCut(paircut); + + + + AliHBTCorrectQInvCorrelFctn* correctqinvCF = new AliHBTCorrectQInvCorrelFctn(); + + analysis->AddTrackFunction(correctqinvCF); + + TFile* outf = TFile::Open("hbtcorrected.root","recreate"); + if (outf == 0x0) + { + cout<<"\n\nERROR: can not open file"<SetInitialValues(0.8,12.0); + correctqinvCF->SetRadiusConvergenceTreshold(0.005); + correctqinvCF->SetLambdaConvergenceTreshold(0.01); + while (1) + { + iter++; + dirname = ietration+iter; + TDirectory* dir = outf->mkdir(dirname,"results after "+dirname); + outf->Write(); + if (dir) dir->cd(); + else + { + cout<<"\n\nERROR: can not make an directory in file named"<Process("Tracks"); + +/*****************************************/ + dir->cd(); + correctqinvCF->WriteAll(); +/*****************************************/ + if (correctqinvCF->IsConverged()) break; + else + { + correctqinvCF->SetInitialValues(correctqinvCF->GetFittedLambda(),correctqinvCF->GetFittedRadius()); + } + } + outf->cd(); + analysis->Write(); + delete outf; +} diff --git a/HBTAN/libHBTAN.pkg b/HBTAN/libHBTAN.pkg index 219216a50d0..6dda49eb980 100644 --- a/HBTAN/libHBTAN.pkg +++ b/HBTAN/libHBTAN.pkg @@ -13,7 +13,7 @@ AliHBTMonDistributionFctns.cxx AliHBTMonResolutionFctns.cxx \ AliHBTLLWeights.cxx AliHBTLLWeightFctn.cxx \ AliHBTLLWeightsPID.cxx AliHBTLLWeightTheorFctn.cxx \ AliHBTPositionRandomizer.cxx AliHBTEventBuffer.cxx \ -AliHBTCorrFitFctn.cxx +AliHBTCorrFitFctn.cxx AliHBTCorrectCorrelFctn.cxx FSRCS = fsiini.F fsiw.F led_bldata.F ltran12.F -- 2.39.3