StRoot: StHbtMaker/CorrFctn/BPLCMSFrame3DCorrFctn

00001 /***************************************************************************
00002  *
00003  * $Id: BPLCMSFrame3DCorrFctn_SIM.cxx,v 1.6 2003/09/02 17:58:20 perev 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  *   THIS IS A SIMULATION CORRELATION FUNCTION CLASS !!!
00011  *
00012  *
00013  ***************************************************************************
00014  *
00015  * $Log: BPLCMSFrame3DCorrFctn_SIM.cxx,v $
00016  * Revision 1.6  2003/09/02 17:58:20  perev
00017  * gcc 3.2 updates + WarnOff
00018  *
00019  * Revision 1.5  2003/01/31 19:21:09  magestro
00020  * Cleared up simple compiler warnings on i386_linux24
00021  *
00022  * Revision 1.4  2001/05/23 00:19:05  lisa
00023  * Add in Smearing classes and methods needed for momentum resolution studies and correction
00024  *
00025  * Revision 1.3  2000/10/08 17:11:07  lisa
00026  * now toggle between Num and Den at BEGINNING of BPLCMSFrame3DCorrFctn_SIM::AddMixedPair()
00027  *
00028  * Revision 1.2  2000/10/05 23:08:59  lisa
00029  * Added kT-dependent radii to mixed-event simulator AND implemented AverageSeparation Cut and CorrFctn
00030  *
00031  * Revision 1.1  2000/09/14 18:36:53  lisa
00032  * Added Qinv and ExitSep pair cuts and BPLCMSFrame3DCorrFctn_SIM CorrFctn
00033  *
00034  *
00035  *
00036  **************************************************************************/
00037 
00038 
00039 #include "StHbtMaker/CorrFctn/BPLCMSFrame3DCorrFctn_SIM.h"
00040 #include "StHbtMaker/Infrastructure/StHbtSmearPair.h"
00041 
00042 #include <cstdio>
00043 
00044 #ifdef __ROOT__ 
00045 ClassImp(BPLCMSFrame3DCorrFctn_SIM)
00046 #endif
00047 
00048 //____________________________
00049 BPLCMSFrame3DCorrFctn_SIM::BPLCMSFrame3DCorrFctn_SIM(char* title, const int& nbins, const float& QLo, const float& QHi){
00050 
00051   // set some stuff...
00052   mQinvNormLo = 0.15;
00053   mQinvNormHi = 0.18;
00054   mNumRealsNorm = 0;
00055   mNumMixedNorm = 0;
00056   mCorrection = 0;  // pointer to Coulomb Correction object
00057 
00058   mPairCut = 0; // added Sept2000 - CorrFctn-specific PairCut
00059 
00060   mToggleNumDen = 0; // this toggles between a pair being given to numerator or denominator.
00061 
00062   mLambda=0.;
00063   mRside2=25.0;
00064   mRout2=25.0;
00065   mRlong2=25.0;
00066 
00067   mRout_alpha = mRside_alpha = mRlong_alpha = 0;  // set these zero by default
00068 
00069   // set up numerator
00070   char TitNum[100] = "Num";
00071   strcat(TitNum,title);
00072   mNumerator = new StHbt3DHisto(TitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00073   // set up denominator
00074   char TitDen[100] = "Den";
00075   strcat(TitDen,title);
00076   mDenominator = new StHbt3DHisto(TitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00077   // set up ratio
00078   char TitRat[100] = "Rat";
00079   strcat(TitRat,title);
00080   mRatio = new StHbt3DHisto(TitRat,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
00081 
00082   // to enable error bar calculation...
00083   mNumerator->Sumw2();
00084   mDenominator->Sumw2();
00085   mRatio->Sumw2();
00086 
00087   // 20feb2001 - add some more histograms that are filled if there is any
00088   //             momentum smearing
00089   
00090   char titInv[100] = "QinvRes";
00091   strcat(titInv,title);
00092   mResolutionHistos[0] = new StHbt2DHisto(titInv,title,10,0.0,0.1,50,-0.025,0.025);
00093   //
00094   char titOut[100] = "QoutRes";
00095   strcat(titOut,title);
00096   mResolutionHistos[1] = new StHbt2DHisto(titOut,title,10,0.0,0.1,50,-0.025,0.025);
00097   //
00098   char titSid[100] = "QsidRes";
00099   strcat(titSid,title);
00100   mResolutionHistos[2] = new StHbt2DHisto(titSid,title,10,0.0,0.1,50,-0.025,0.025);
00101   //
00102   char titLon[100] = "QlonRes";
00103   strcat(titLon,title);
00104   mResolutionHistos[3] = new StHbt2DHisto(titLon,title,10,0.0,0.1,50,-0.025,0.025);
00105 
00106 }
00107 
00108 //____________________________
00109 BPLCMSFrame3DCorrFctn_SIM::~BPLCMSFrame3DCorrFctn_SIM(){
00110   delete mNumerator;
00111   delete mDenominator;
00112   delete mRatio;
00113 }
00114 //_________________________
00115 void BPLCMSFrame3DCorrFctn_SIM::Finish(){
00116 
00117   mRatio->Divide(mNumerator,mDenominator);
00118 }
00119 
00120 //____________________________
00121 StHbtString BPLCMSFrame3DCorrFctn_SIM::Report(){
00122   string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n";
00123   char ctemp[100];
00124   sprintf(ctemp,"    THIS IS A SIMULATION CORRELATION FUNCTION CLASS!!!\n");
00125   stemp += ctemp;
00126   sprintf(ctemp,"Number of entries in numerator:\t%E\n",mNumerator->GetEntries());
00127   stemp += ctemp;
00128   sprintf(ctemp,"Number of entries in denominator:\t%E\n",mDenominator->GetEntries());
00129   stemp += ctemp;
00130   sprintf(ctemp,"Number of entries in ratio:\t%E\n",mRatio->GetEntries());
00131   stemp += ctemp;
00132   sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",mQinvNormLo,mQinvNormHi);
00133   stemp += ctemp;
00134   sprintf(ctemp,"Number of pairs in Normalization region was:\n");
00135   stemp += ctemp;
00136   sprintf(ctemp,"In numerator:\t%lu\t In denominator:\t%lu\n",mNumRealsNorm,mNumMixedNorm);
00137   stemp += ctemp;
00138   if (mCorrection)
00139     {
00140       float radius = mCorrection->GetRadius();
00141       sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius);
00142     }
00143   else
00144     {
00145       sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n");
00146     }
00147   stemp += ctemp;
00148 
00149   if (mPairCut){
00150     sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
00151     stemp += ctemp;
00152     stemp += mPairCut->Report();
00153   }
00154   else{
00155     sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
00156     stemp += ctemp;
00157   }
00158 
00159   //  
00160   StHbtString returnThis = stemp;
00161   return returnThis;
00162 }
00163 //____________________________
00164 void BPLCMSFrame3DCorrFctn_SIM::AddRealPair(const StHbtPair* pair){
00165 
00166   // note that in this SIMULATION correlation function class, we do NOTHING with real
00167   // pairs.  The numerator is just a weighted denominator, filled in AddMixedPair (below)
00168 
00169   /* no-op */
00170   return;
00171 
00172 }
00173 
00174 
00175 
00176 
00177 //____________________________
00178 void BPLCMSFrame3DCorrFctn_SIM::AddMixedPair(const StHbtPair* pair){
00179 
00180   mToggleNumDen = !mToggleNumDen;  // toggle at the BEGINNING of call instead of end
00181                                    // so that SAME pairs go into Num and Den as long as they pass cut
00182 
00183 
00184   if (mPairCut){
00185     if (!(mPairCut->Pass(pair))){
00186 //       // mike is checking stuff here
00187 //       double qOut = fabs(pair->qOutCMS());
00188 //       double qSide = fabs(pair->qSideCMS());
00189 //       double qLong = fabs(pair->qLongCMS());
00190 //       //      if ((qSide > 0.035)||(qLong > 0.035)){
00191 //       if ((qSide > 0.0)||(qLong > 0.0)){
00192 //      cout << qOut << " " << qSide << " " << qLong << " " << 
00193 //        pair->NominalTpcEntranceSeparation() << " " <<
00194 //        pair->track1()->NominalTpcEntrancePoint() << " " <<
00195 //        pair->track2()->NominalTpcEntrancePoint() << endl;
00196 //      int ijunk;
00197 //      cout << "Enter a junk integer : ";
00198 //      cin >> ijunk;
00199 //       }
00200 //      // end of checking stuff
00201       return;
00202     }
00203   }
00204 
00205   double weight=1.0;
00206 
00207   // these are the values we use to assign correlation WEIGHT
00208   double qOut = fabs(pair->qOutCMS());
00209   double qSide = fabs(pair->qSideCMS());
00210   double qLong = fabs(pair->qLongCMS());
00211 
00212   // these are the values we BIN in - same as above if no smearing turned on...
00213   double qOut_bin,qSide_bin,qLong_bin;  // scope 'em!
00214   if (mSmearPair){
00215     mSmearPair->SetUnsmearedPair(pair);
00216     qOut_bin = fabs(mSmearPair->SmearedPair().qOutCMS());
00217     qSide_bin = fabs(mSmearPair->SmearedPair().qSideCMS());
00218     qLong_bin = fabs(mSmearPair->SmearedPair().qLongCMS());
00219     // fill resolution histos...
00220     double Qinv = fabs(pair->qInv());
00221     if (Qinv < 0.1){
00222       double Qinv_smear = fabs(mSmearPair->SmearedPair().qInv());
00223       mResolutionHistos[0]->Fill(Qinv,Qinv_smear-Qinv);
00224       mResolutionHistos[1]->Fill(qOut,qOut_bin-qOut);
00225       mResolutionHistos[2]->Fill(qSide,qSide_bin-qSide);
00226       mResolutionHistos[3]->Fill(qLong,qLong_bin-qLong);
00227     }
00228   }
00229   else{
00230     qOut_bin  = qOut;
00231     qSide_bin = qSide;
00232     qLong_bin = qLong;
00233   }
00234 
00235 
00236   // note the Coulomb bit below...
00237   // for the *suppression* ("numerator" pairs), we use the TRUE momentum, as Nature would
00238   // for the *correction* ("denominator" pairs), we use the SMEARED momentum, as an experimenter would
00239 
00240   if (mToggleNumDen){
00241     if (mCorrection){
00242       if (mSmearPair){
00243         weight = mCorrection->CoulombCorrect(&(mSmearPair->SmearedPair()));
00244       }
00245       else{
00246         weight = mCorrection->CoulombCorrect(pair);
00247       }
00248     }
00249     mDenominator->Fill(qOut_bin,qSide_bin,qLong_bin,weight);
00250   }
00251   else{
00252 
00253     if (mCorrection) weight = mCorrection->CoulombCorrect(pair);
00254 
00255     // update 1oct2000 - now allow mT dependent evolution of Radii
00256 
00257     float mT2 = (pair->fourMomentumSum().mt2())/4.0;  // divide by 4.0 because want mT of (p1+p2)/2 
00258                                                        // and factor of two gets squared in mt2
00259 
00260     float Rout2,Rside2,Rlong2;
00261 
00262     if (mRout_alpha!=0){Rout2 = mRout2*::pow(mT2,mRout_alpha);}
00263     else{Rout2 = mRout2;}
00264 
00265     if (mRside_alpha!=0){Rside2 = mRside2*::pow(mT2,mRside_alpha);}
00266     else{Rside2 = mRside2;}
00267 
00268     if (mRlong_alpha!=0){Rlong2 = mRlong2*::pow(mT2,mRlong_alpha);}
00269     else{Rlong2 = mRlong2;}
00270       
00271     double CorrWeight = 1.0 + 
00272       mLambda*exp((-qOut*qOut*Rout2 -qSide*qSide*Rside2 -qLong*qLong*Rlong2)/0.038936366329);
00273     CorrWeight *= weight;
00274     mNumerator->Fill(qOut_bin,qSide_bin,qLong_bin,CorrWeight);
00275   }
00276 
00277 }
00278 
00279 
00280
StHbtMaker/CorrFctn/BPLCMSFrame3DCorrFctn_SIM.cxx
