StRoot: StHbtMaker/CorrFctn/BPLCMSFrame3DCorrFctn.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: BPLCMSFrame3DCorrFctn.cxx,v 1.7 2003/01/31 19:20:54 magestro Exp $
00004  *
00005  * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
00006  ***************************************************************************
00007  *
00008  * Description: part of STAR HBT Framework: StHbtMaker package
00009  *   3D Bertsch-Pratt decomposition in the LCMS frame
00010  *
00011  ***************************************************************************
00012  *
00013  * $Log: BPLCMSFrame3DCorrFctn.cxx,v $
00014  * Revision 1.7  2003/01/31 19:20:54  magestro
00015  * Cleared up simple compiler warnings on i386_linux24
00016  *
00017  * Revision 1.6  2002/06/07 22:51:38  lisa
00018  * Widely used BPLCMSFrame3DCorrFctn class now accumulates UNcorrected denominator and has a WriteOutHistos method
00019  *
00020  * Revision 1.5  2001/05/23 00:19:04  lisa
00021  * Add in Smearing classes and methods needed for momentum resolution studies and correction
00022  *
00023  * Revision 1.4  2000/10/26 19:48:49  rcwells
00024  * Added functionality for Coulomb correction of <qInv> in 3D correltions
00025  *
00026  * Revision 1.3  2000/09/14 18:36:53  lisa
00027  * Added Qinv and ExitSep pair cuts and BPLCMSFrame3DCorrFctn_SIM CorrFctn
00028  *
00029  * Revision 1.2  2000/08/23 19:43:43  lisa
00030  * added alternate normalization algorithm to 3d CorrFctns in case normal one fails
00031  *
00032  * Revision 1.1  2000/08/17 20:48:39  lisa
00033  * Adding correlationfunction in LCMS frame
00034  *
00035  *
00036  **************************************************************************/
00037 
00038 #include "StHbtMaker/CorrFctn/BPLCMSFrame3DCorrFctn.h"
00039 //#include "StHbtMaker/Infrastructure/StHbtHisto.hh"
00040 #include <cstdio>
00041 
00042 #ifdef __ROOT__ 
00043 ClassImp(BPLCMSFrame3DCorrFctn)
00044 #endif
00045 
00046 //____________________________
00047 BPLCMSFrame3DCorrFctn::BPLCMSFrame3DCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi){
00048 
00049   // set some stuff...
00050   mQinvNormLo = 0.15;
00051   mQinvNormHi = 0.18;
00052   mNumRealsNorm = 0;
00053   mNumMixedNorm = 0;
00054   mCorrection = 0;  // pointer to Coulomb Correction object
00055 
00056   mPairCut = 0; // added Sept2000 - CorrFctn-specific PairCut
00057 
00058   mSmearPair = 0; // no resolution correction unless user sets SmearPair
00059 
00060   // set up numerator
00061   char TitNum[100] = "Num";
00062   strcat(TitNum,title);
00063   mNumerator = new StHbt3DHisto(TitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00064   // set up denominator
00065   char TitDen[100] = "Den";
00066   strcat(TitDen,title);
00067   mDenominator = new StHbt3DHisto(TitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00068   // set up uncorrected denominator
00069   char TitDenUncoul[100] = "DenNoCoul";
00070   strcat(TitDenUncoul,title);
00071   mUncorrectedDenominator = new StHbt3DHisto(TitDenUncoul,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00072   // set up ratio
00073   char TitRat[100] = "Rat";
00074   strcat(TitRat,title);
00075   mRatio = new StHbt3DHisto(TitRat,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00076   // set up ave qInv
00077   char TitQinv[100] = "Qinv";
00078   strcat(TitQinv,title);
00079   mQinvHisto = new StHbt3DHisto(TitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00080 
00081   // to enable error bar calculation...
00082   mNumerator->Sumw2();
00083   mDenominator->Sumw2();
00084   mUncorrectedDenominator->Sumw2();
00085   mRatio->Sumw2();
00086 
00087   // Following histos are for the momentum resolution correction
00088   // they are filled only if a StHbtSmear object is plugged in
00089   // here comes the "idea" numerator and denominator and ratio...
00090   char TitNumID[100] = "IDNum";
00091   strcat(TitNumID,title);
00092   mIDNumHisto = new StHbt3DHisto(TitNumID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00093   char TitDenID[100] = "IDDen";
00094   strcat(TitDenID,title);
00095   mIDDenHisto = new StHbt3DHisto(TitDenID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00096   char TitRatID[100] = "IDRat";
00097   strcat(TitRatID,title);
00098   mIDRatHisto = new StHbt3DHisto(TitRatID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00099 
00100   mIDNumHisto->Sumw2();
00101   mIDDenHisto->Sumw2();
00102   mIDRatHisto->Sumw2();
00103 
00104   //
00105   // here comes the "smeared" numerator and denominator...
00106   char TitNumSM[100] = "SMNum";
00107   strcat(TitNumSM,title);
00108   mSMNumHisto = new StHbt3DHisto(TitNumSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00109   char TitDenSM[100] = "SMDen";
00110   strcat(TitDenSM,title);
00111   mSMDenHisto = new StHbt3DHisto(TitDenSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00112   char TitRatSM[100] = "SMRat";
00113   strcat(TitRatSM,title);
00114   mSMRatHisto = new StHbt3DHisto(TitRatSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00115   //
00116   mSMNumHisto->Sumw2();
00117   mSMDenHisto->Sumw2();
00118   mSMRatHisto->Sumw2();
00119   //
00120   // here comes the correction factor (which is just ratio of ideal ratio to smeared ratio)
00121   char TitCorrection[100] = "CorrectionFactor";
00122   strcat(TitCorrection,title);
00123   mCorrectionHisto = new StHbt3DHisto(TitCorrection,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);  
00124   mCorrectionHisto->Sumw2();
00125   // here comes the fully corrected correlation function
00126   char TitCorrCF[100] = "CorrectedCF";
00127   strcat(TitCorrCF,title);
00128   mCorrCFHisto = new StHbt3DHisto(TitCorrCF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00129   mCorrCFHisto->Sumw2();
00130 
00131   // user can (and should) override these defaults...
00132   mLambda = 0.6;
00133   mRout2 = 6.0*6.0;
00134   mRside2 = 6.0*6.0;
00135   mRlong2 = 7.0*7.0;
00136 
00137 }
00138 
00139 //____________________________
00140 BPLCMSFrame3DCorrFctn::~BPLCMSFrame3DCorrFctn(){
00141   delete mNumerator;
00142   delete mDenominator;
00143   delete mRatio;
00144   delete mQinvHisto;
00145   delete mIDNumHisto;
00146   delete mIDDenHisto;
00147   delete mIDRatHisto;
00148   delete mSMNumHisto;
00149   delete mSMDenHisto;
00150   delete mSMRatHisto;
00151   delete mCorrectionHisto;
00152   delete mCorrCFHisto;
00153 }
00154 
00155 //_________________________
00156 void BPLCMSFrame3DCorrFctn::WriteOutHistos(){
00157 
00158   mNumerator->Write();
00159   mDenominator->Write();
00160   mUncorrectedDenominator->Write();
00161   mRatio->Write();
00162   mQinvHisto->Write();
00163 
00164   if (mSmearPair){
00165     mIDNumHisto->Write();
00166     mIDDenHisto->Write();
00167     mIDRatHisto->Write();
00168     //
00169     mSMNumHisto->Write();
00170     mSMDenHisto->Write();
00171     mSMRatHisto->Write();
00172     //
00173     mCorrectionHisto->Write();
00174     mCorrCFHisto->Write();
00175   }
00176 }
00177 
00178 //_________________________
00179 void BPLCMSFrame3DCorrFctn::Finish(){
00180   // here is where we should normalize, fit, etc...
00181   double NumFact,DenFact;
00182   if ((mNumRealsNorm !=0) && (mNumMixedNorm !=0)){
00183     NumFact = double(mNumRealsNorm);
00184     DenFact = double(mNumMixedNorm);
00185   }
00186   // can happen that the mNumRealsNorm and mNumMixedNorm = 0 if you do non-standard
00187   //   things like making a new CorrFctn and just setting the Numerator and Denominator
00188   //   from OTHER CorrFctns which you read in (like when doing parallel processing) 
00189   else{
00190     cout << "Warning! - no normalization constants defined - I do the best I can..." << endl;
00191     int nbins = mNumerator->GetNbinsX();
00192     int half_way = nbins/2;
00193     NumFact = mNumerator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
00194     DenFact = mDenominator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
00195   }
00196 
00197   mRatio->Divide(mNumerator,mDenominator,DenFact,NumFact);
00198   mQinvHisto->Divide(mUncorrectedDenominator);
00199 
00200   // now do all the resolution correction stuff..
00201   if (mSmearPair){  // but only do it if we have been working with a SmearPair
00202     mIDRatHisto->Divide(mIDNumHisto,mIDDenHisto);
00203     mSMRatHisto->Divide(mSMNumHisto,mSMDenHisto);
00204     mCorrectionHisto->Divide(mIDRatHisto,mSMRatHisto);
00205     mCorrCFHisto->Multiply(mRatio,mCorrectionHisto);
00206   }
00207 
00208 }
00209 
00210 //____________________________
00211 StHbtString BPLCMSFrame3DCorrFctn::Report(){
00212   string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n";
00213   char ctemp[100];
00214   sprintf(ctemp,"Number of entries in numerator:\t%E\n",mNumerator->GetEntries());
00215   stemp += ctemp;
00216   sprintf(ctemp,"Number of entries in denominator:\t%E\n",mDenominator->GetEntries());
00217   stemp += ctemp;
00218   sprintf(ctemp,"Number of entries in ratio:\t%E\n",mRatio->GetEntries());
00219   stemp += ctemp;
00220   sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",mQinvNormLo,mQinvNormHi);
00221   stemp += ctemp;
00222   sprintf(ctemp,"Number of pairs in Normalization region was:\n");
00223   stemp += ctemp;
00224   sprintf(ctemp,"In numerator:\t%lu\t In denominator:\t%lu\n",mNumRealsNorm,mNumMixedNorm);
00225   stemp += ctemp;
00226   if (mCorrection)
00227     {
00228       float radius = mCorrection->GetRadius();
00229       sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius);
00230     }
00231   else
00232     {
00233       sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n");
00234     }
00235   stemp += ctemp;
00236 
00237   if (mPairCut){
00238     sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
00239     stemp += ctemp;
00240     stemp += mPairCut->Report();
00241   }
00242   else{
00243     sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
00244     stemp += ctemp;
00245   }
00246 
00247   //  
00248   StHbtString returnThis = stemp;
00249   return returnThis;
00250 }
00251 //____________________________
00252 void BPLCMSFrame3DCorrFctn::AddRealPair(const StHbtPair* pair){
00253 
00254   if (mPairCut){
00255     if (!(mPairCut->Pass(pair))) return;
00256   }
00257 
00258   double Qinv = fabs(pair->qInv());   // note - qInv() will be negative for identical pairs...
00259   if ((Qinv < mQinvNormHi) && (Qinv > mQinvNormLo)) mNumRealsNorm++;
00260   double qOut = fabs(pair->qOutCMS());
00261   double qSide = fabs(pair->qSideCMS());
00262   double qLong = fabs(pair->qLongCMS());
00263 
00264   mNumerator->Fill(qOut,qSide,qLong);
00265 }
00266 //____________________________
00267 void BPLCMSFrame3DCorrFctn::AddMixedPair(const StHbtPair* pair){
00268 
00269   if (mPairCut){
00270     if (!(mPairCut->Pass(pair))) return;
00271   }
00272 
00273   double CoulombWeight = (mCorrection ? mCorrection->CoulombCorrect(pair) : 1.0);
00274 
00275   double Qinv = fabs(pair->qInv());   // note - qInv() will be negative for identical pairs...
00276   if ((Qinv < mQinvNormHi) && (Qinv > mQinvNormLo)) mNumMixedNorm++;
00277   double qOut = fabs(pair->qOutCMS());
00278   double qSide = fabs(pair->qSideCMS());
00279   double qLong = fabs(pair->qLongCMS());
00280 
00281   mDenominator->Fill(qOut,qSide,qLong,CoulombWeight);
00282   mUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0);
00283   mQinvHisto->Fill(qOut,qSide,qLong,Qinv);
00284 
00285   // now for the momentum resolution stuff...
00286   if (mSmearPair){
00287     double CorrWeight =  1.0 + 
00288       mLambda*exp((-qOut*qOut*mRout2 -qSide*qSide*mRside2 -qLong*qLong*mRlong2)/0.038936366329);
00289     CorrWeight *= CoulombWeight;  // impt.
00290 
00291     mIDNumHisto->Fill(qOut,qSide,qLong,CorrWeight);
00292     mIDDenHisto->Fill(qOut,qSide,qLong,CoulombWeight);
00293 
00294     mSmearPair->SetUnsmearedPair(pair);
00295     double qOut_prime = fabs(mSmearPair->SmearedPair().qOutCMS());
00296     double qSide_prime = fabs(mSmearPair->SmearedPair().qSideCMS());
00297     double qLong_prime = fabs(mSmearPair->SmearedPair().qLongCMS());
00298 
00299     mSMNumHisto->Fill(qOut_prime,qSide_prime,qLong_prime,CorrWeight);
00300 
00301     double SmearedCoulombWeight = ( mCorrection ? 
00302                                     mCorrection->CoulombCorrect(&(mSmearPair->SmearedPair())) : 
00303                                     1.0);
00304 
00305     mSMDenHisto->Fill(qOut_prime,qSide_prime,qLong_prime,SmearedCoulombWeight);
00306   }
00307 }
00308