StRoot: StHbtMaker/CorrFctn/YKPLongitudinal.cxx Source File

00001 /***************************************************************************
00002  *
00003  * Author: Dominik Flierl, flierl@bnl.gov
00004  ***************************************************************************
00005  *
00006  * Description: 
00007  *  YKP parametrization in longitudinal direction in kt-Y-Bins
00008  *
00009  *
00010  **************************************************************************/
00011 
00012 
00013 #include "StHbtMaker/CorrFctn/YKPLongitudinal.h"
00014 #include <cstdio>
00015 
00016 #ifdef __ROOT__ 
00017 ClassImp(YKPLongitudinal)
00018 #endif
00019 
00020 //____________________________
00021 YKPLongitudinal::YKPLongitudinal( TString ctype, TString frame , 
00022                                   const int& nbins, const float& QLo, const float& QHi, 
00023                                   const int& nKtbins, const double& ktMin, const double& ktMax,
00024                                   const int& nYbins,  const double& YMin,  const double& YMax,
00025                                   const int& nbinsQINV, const float& QLoQINV , const float& QHiQINV )
00026 {
00027   // set correlation function type
00028   if ( !ctype.CompareTo("YKP") || !ctype.CompareTo("BP") )
00029     {
00030       mCtype = ctype  ;
00031       cout << "correlation function set to :" << ctype << endl ;
00032     }
00033   else 
00034     {
00035       cout << "Error: correlationfunction didn't fit BP or YKP -> YKP used !\n" ;
00036       mCtype = "YKP" ;
00037     }
00038     
00039   // set frame
00040   if ( !frame.CompareTo("LCMS") || !frame.CompareTo("CMS") || !frame.CompareTo("PRF")) 
00041     {
00042       mFrame = frame  ;
00043       cout << "frame set to :" << frame << endl ;
00044     }
00045   else 
00046     {
00047       cout << "Error: frame choice didn't fit LCMS,CMS or PRF -> LCMS used !\n" ;
00048       mFrame = "LCMS" ;
00049     }
00050   
00051   // get intializers
00052   mNumberKtBins = nKtbins ;
00053   mNumberYBins = nYbins ;
00054     
00055   // note : bin = binKt + (binY-1) * mNumberKtBins
00056   mNumberBins = mNumberKtBins * mNumberYBins ;
00057   
00058   // pointer to Coulomb Correction object must be set from macro
00059   mCorrection = 0 ;  
00060   
00061   // set up 3d histos
00062   mNumerator    = new StHbt3DHisto[mNumberBins] ;
00063   mDenominator  = new StHbt3DHisto[mNumberBins] ;
00064   mRatio        = new StHbt3DHisto[mNumberBins] ;
00065   mQinvNumerator   = new StHbt1DHisto[mNumberBins] ;
00066   mQinvDenominator = new StHbt1DHisto[mNumberBins] ;
00067   mQinvRatio       = new StHbt1DHisto[mNumberBins] ;
00068 
00069   char TitNum[100] ; 
00070   char TitDen[100] ; 
00071   char TitRat[100] ;
00072   char TitQinvNum[100] ;
00073   char TitQinvDen[100] ;
00074   char TitQinvRatio[100] ;
00075 
00076   char xAxisTitle[100] ; 
00077   char yAxisTitle[100] ; 
00078   char zAxisTitle[100] ;
00079   if (mCtype == "YKP")
00080     {
00081       sprintf(xAxisTitle,"qLong") ;
00082       sprintf(yAxisTitle,"qPerp") ;
00083       sprintf(zAxisTitle,"qNull") ;
00084     }
00085   else if (mCtype == "BP")
00086     {
00087       sprintf(xAxisTitle,"qSide") ;
00088       sprintf(yAxisTitle,"qOut") ;
00089       sprintf(zAxisTitle,"qLong") ;
00090     }
00091 
00092   for(int ktindex = 0 ; ktindex < mNumberKtBins ; ktindex++)
00093     { 
00094       for(int yindex = 0 ; yindex < mNumberYBins ; yindex++)    
00095         { 
00096           // numerator
00097           sprintf(TitNum,"NumKt%dY%d",ktindex+1,yindex+1) ;
00098           mNumerator[ktindex+yindex*mNumberKtBins].SetName(TitNum);
00099           mNumerator[ktindex+yindex*mNumberKtBins].SetTitle(TitNum);
00100           mNumerator[ktindex+yindex*mNumberKtBins].SetBins(nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi) ;
00101           mNumerator[ktindex+yindex*mNumberKtBins].SetXTitle(xAxisTitle); 
00102           mNumerator[ktindex+yindex*mNumberKtBins].SetYTitle(yAxisTitle);  
00103           mNumerator[ktindex+yindex*mNumberKtBins].SetZTitle(zAxisTitle);
00104           // denominator
00105           sprintf(TitDen,"DenKt%dY%d",ktindex+1,yindex+1) ;
00106           mDenominator[ktindex+yindex*mNumberKtBins].SetName(TitDen);
00107           mDenominator[ktindex+yindex*mNumberKtBins].SetTitle(TitDen);
00108           mDenominator[ktindex+yindex*mNumberKtBins].SetBins(nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi) ;
00109           mDenominator[ktindex+yindex*mNumberKtBins].SetXTitle(xAxisTitle); 
00110           mDenominator[ktindex+yindex*mNumberKtBins].SetYTitle(yAxisTitle);  
00111           mDenominator[ktindex+yindex*mNumberKtBins].SetZTitle(zAxisTitle);
00112           // ratio
00113           sprintf(TitRat,"RatKt%dY%d",ktindex+1,yindex+1) ;
00114           mRatio[ktindex+yindex*mNumberKtBins].SetName(TitRat);
00115           mRatio[ktindex+yindex*mNumberKtBins].SetTitle(TitRat);
00116           mRatio[ktindex+yindex*mNumberKtBins].SetBins(nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi) ;
00117           mRatio[ktindex+yindex*mNumberKtBins].SetXTitle(xAxisTitle); 
00118           mRatio[ktindex+yindex*mNumberKtBins].SetYTitle(yAxisTitle);  
00119           mRatio[ktindex+yindex*mNumberKtBins].SetZTitle(zAxisTitle);
00120           // qinv numerator
00121           sprintf(TitQinvNum,"QinvNumeratorKt%dY%d",ktindex+1,yindex+1) ;
00122           mQinvNumerator[ktindex+yindex*mNumberKtBins].SetName(TitQinvNum);
00123           mQinvNumerator[ktindex+yindex*mNumberKtBins].SetTitle(TitQinvNum);
00124           mQinvNumerator[ktindex+yindex*mNumberKtBins].SetBins(nbinsQINV,QLoQINV,QHiQINV) ;
00125           mQinvNumerator[ktindex+yindex*mNumberKtBins].SetXTitle("qInv"); 
00126           // qinv denominator
00127           sprintf(TitQinvDen,"QinvDenominatorKt%dY%d",ktindex+1,yindex+1) ;
00128           mQinvDenominator[ktindex+yindex*mNumberKtBins].SetName(TitQinvDen);
00129           mQinvDenominator[ktindex+yindex*mNumberKtBins].SetTitle(TitQinvDen);
00130           mQinvDenominator[ktindex+yindex*mNumberKtBins].SetBins(nbinsQINV,QLoQINV,QHiQINV) ;
00131           mQinvDenominator[ktindex+yindex*mNumberKtBins].SetXTitle("qInv"); 
00132           // qinv ratio
00133           sprintf(TitQinvRatio,"QinvRatioKt%dY%d",ktindex+1,yindex+1) ;
00134           mQinvRatio[ktindex+yindex*mNumberKtBins].SetName(TitQinvRatio);
00135           mQinvRatio[ktindex+yindex*mNumberKtBins].SetTitle(TitQinvRatio);
00136           mQinvRatio[ktindex+yindex*mNumberKtBins].SetBins(nbinsQINV,QLoQINV,QHiQINV) ;
00137           mQinvRatio[ktindex+yindex*mNumberKtBins].SetXTitle("qInv") ; 
00138           // to enable error bar calculation...
00139           mNumerator[ktindex+yindex*mNumberKtBins].Sumw2() ;
00140           mDenominator[ktindex+yindex*mNumberKtBins].Sumw2() ;
00141           mRatio[ktindex+yindex*mNumberKtBins].Sumw2() ;
00142           mQinvNumerator[ktindex+yindex*mNumberKtBins].Sumw2() ;
00143           mQinvDenominator[ktindex+yindex*mNumberKtBins].Sumw2() ; 
00144           mQinvRatio[ktindex+yindex*mNumberKtBins].Sumw2() ;
00145         }
00146     }
00147 
00148 
00150   // set up y pt bins boundaries
00152   // kt
00153   mktBinsMin = new double[mNumberKtBins] ;
00154   mktBinsMax = new double[mNumberKtBins] ;
00155   double ktstep = (double)(ktMax-ktMin)/(double) mNumberKtBins ;
00156   for(int ktindex = 0 ; ktindex < mNumberKtBins  ;  ktindex++ )
00157     {
00158       mktBinsMin[ktindex] = ((double) ktindex) * ktstep + ktMin ;
00159       mktBinsMax[ktindex] = ((double) (ktindex+1)) * ktstep + ktMin ;
00160     }
00161   // rap
00162   mYBinsMin  = new double[mNumberKtBins] ;
00163   mYBinsMax  = new double[mNumberYBins] ;
00164   double ystep  = (double)(YMax-YMin)/(double) mNumberYBins ;
00165   for(int yindex = 0 ; yindex < mNumberYBins  ;  yindex++ )
00166     {
00167       mYBinsMin[yindex] = ((double) yindex) * ystep + YMin ;
00168       mYBinsMax[yindex] = ((double) (yindex+1)) * ystep + YMin ;
00169     }
00170   
00172   // set up normalization
00174   // set range normalization
00175   mQinvNormLo = 0.15 ;
00176   mQinvNormHi = 0.18 ;
00177   // set up  array with normalization values    
00178   mNumRealsNorm = new unsigned long int [mNumberBins] ;
00179   mNumMixedNorm = new unsigned long int [mNumberBins] ;
00180   for( int index = 0 ; index < mNumberBins ; index++ ) 
00181     {
00182       mNumRealsNorm[index] = mNumMixedNorm[index] = 0 ;
00183     }
00184 }
00185 //____________________________
00186 YKPLongitudinal::~YKPLongitudinal()
00187 {
00188   delete [] mNumRealsNorm ;
00189   delete [] mNumMixedNorm ;
00190   delete [] mNumerator ;
00191   delete [] mDenominator ;
00192   delete [] mRatio ; 
00193   delete [] mQinvNumerator ; 
00194   delete [] mQinvDenominator ; 
00195   delete [] mQinvRatio ;
00196 }
00197 //_________________________
00198 void YKPLongitudinal::Finish()
00199 {
00200   // here is where we can normalize, fit, etc...
00201   double NumFact,DenFact;
00202   for(int ktindex = 0 ; ktindex < mNumberKtBins ; ktindex++)
00203     { 
00204       for(int yindex = 0 ; yindex < mNumberYBins ; yindex++)
00205         { 
00206           // get the normalisation 
00207           if ((mNumRealsNorm[ktindex+yindex*mNumberKtBins] !=0) && (mNumMixedNorm[ktindex+yindex*mNumberKtBins] !=0))
00208             {
00209               NumFact = double(mNumRealsNorm[ktindex+yindex*mNumberKtBins]) ;
00210               DenFact = double(mNumMixedNorm[ktindex+yindex*mNumberKtBins]) ;
00211               cout << " Normalizing with Num/Denom = norm : " << NumFact << "\t" << DenFact << "\t" <<  NumFact/DenFact 
00212                    << " for kt bin : " << ktindex << "   y bin : " << yindex << endl ; 
00213             }
00214           // error no normalization filled
00215           else
00216             {
00217               cout << "Warning! - no normalization possible ..." << endl;
00218               NumFact = 1.0 ;
00219               DenFact = 10.0 ;
00220             }
00221           // normalize ratio
00222           mRatio[ktindex+yindex*mNumberKtBins].Divide(&mNumerator[ktindex+yindex*mNumberKtBins],
00223                                                       &mDenominator[ktindex+yindex*mNumberKtBins],DenFact,NumFact) ;
00224           
00225           // produce qinv ratio
00226           mQinvRatio[ktindex+yindex*mNumberKtBins].Divide(&mQinvNumerator[ktindex+yindex*mNumberKtBins],
00227                                                           &mQinvDenominator[ktindex+yindex*mNumberKtBins]) ;
00228         }
00229     }
00230 }
00231 //____________________________
00232 StHbtString YKPLongitudinal::Report()
00233 {
00234   // report called by manager per standard
00235   // get the overall values
00236   int mNumeratorEntriesAll = 0 ;
00237   int mDenominatorEntriesAll = 0 ;
00238   int mRatioEntriesAll = 0 ;
00239   int mNumRealsNormAll = 0 ;
00240   int mNumMixedNormAll = 0 ;
00241   
00242   for(int ktindex = 0 ; ktindex < mNumberKtBins ; ktindex++)
00243     { 
00244       for(int yindex = 0 ; yindex < mNumberYBins ; yindex++)
00245         { 
00246           mNumeratorEntriesAll   += (int)mNumerator[ktindex+yindex*mNumberKtBins].GetEntries() ;
00247           mDenominatorEntriesAll += (int)mDenominator[ktindex+yindex*mNumberKtBins].GetEntries() ;
00248           mRatioEntriesAll       += (int)mRatio[ktindex+yindex*mNumberKtBins].GetEntries() ;    
00249           mNumRealsNormAll  += mNumRealsNorm[ktindex+yindex*mNumberKtBins] ;
00250           mNumMixedNormAll  += mNumMixedNorm[ktindex+yindex*mNumberKtBins] ;
00251         }
00252     }
00253   // put the report together
00254   string stemp = "YKP Function Report:\n";
00255   char ctemp[100];
00256   sprintf(ctemp,"Number of entries in numerator in all %d histos:\t%d\n",mNumberBins,mNumeratorEntriesAll);
00257   stemp += ctemp;
00258   sprintf(ctemp,"Number of entries in denominator in all %d histos:\t%d\n",mNumberBins,mDenominatorEntriesAll);
00259   stemp += ctemp;
00260   sprintf(ctemp,"Number of entries in ratio in all %d histos:\t%d\n",mNumberBins,mRatioEntriesAll);
00261   stemp += ctemp;
00262   sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",mQinvNormLo,mQinvNormHi);
00263   stemp += ctemp;
00264   sprintf(ctemp,"Number of pairs in Normalization region was:\n");
00265   stemp += ctemp;
00266   sprintf(ctemp,"In numerator:\t%u\t In denominator:\t%u\n",mNumRealsNormAll,mNumMixedNormAll);
00267   stemp += ctemp ;
00268   if (mCorrection)
00269     {
00270       float radius = mCorrection->GetRadius();
00271           sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius) ;
00272     }
00273   else
00274     {
00275       sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n") ;
00276     }
00277   stemp += ctemp ;
00278   
00279   //  return it
00280   StHbtString returnThis = stemp ;
00281   return returnThis ;
00282 }
00283 //____________________________
00284 void YKPLongitudinal::AddRealPair(const StHbtPair* pair)
00285 {
00286   // get the pair momenta 
00287   double mKt = pair->kT() ; 
00288   double mRap = pair->rap() ;
00289   double mQinv = fabs(pair->qInv()) ; 
00290   // sort it into the approriate histo
00291   for(int ktindex = 0 ; ktindex < mNumberKtBins ; ktindex++)
00292     { 
00293       if( mKt >= mktBinsMin[ktindex] && mKt < mktBinsMax[ktindex] )
00294         {
00295           for(int yindex = 0 ; yindex < mNumberYBins ; yindex++)
00296             { 
00297               if( mRap >= mYBinsMin[yindex] && mRap < mYBinsMax[yindex] ) 
00298                 {
00299                   // sum up qinv pairs to do later the normalisation 
00300                   if ((mQinv < mQinvNormHi) && (mQinv > mQinvNormLo)) mNumRealsNorm[ktindex+yindex*mNumberKtBins]++;
00301                   // fill 1d qinv histo
00302                   mQinvNumerator[ktindex+yindex*mNumberKtBins].Fill(mQinv) ;
00303                   // fill 3d : choose correlationfunction type and  frame 
00304                   if (mCtype == "YKP")
00305                     {
00306                       double qlong, qperp, q0 ;
00307                       if(mFrame =="LCMS") { pair->qYKPLCMS(qlong,qperp,q0) ;}
00308                       else if(mFrame =="CMS") { pair->qYKPCMS(qlong,qperp,q0) ;}
00309                       else if(mFrame =="PRF") { pair->qYKPPF(qlong,qperp,q0) ;}
00310                       // fill 3d histo
00311                       mNumerator[ktindex+yindex*mNumberKtBins].Fill(qlong,qperp,q0) ;
00312                     }
00313                   else if (mCtype == "BP")
00314                     {
00315                       double qside, qout, qlong;
00316                       if(mFrame =="LCMS")     { qside = pair->qSideCMS(); qout = pair->qOutCMS(); qlong = pair->qLongCMS(); }
00317                       else if(mFrame =="CMS") { qside = pair->qSideBf();  qout = pair->qOutBf();  qlong = pair->qLongBf();  }
00318                       else if(mFrame =="PRF") { qside = pair->qSidePf();  qout = pair->qOutPf();  qlong = pair->qLongPf(); }
00319                       // fill 3d histo
00320                       mNumerator[ktindex+yindex*mNumberKtBins].Fill(fabs(qside),fabs(qout),fabs(qlong)) ;
00321                     }
00322                   // end yindex loop
00323                   break ; 
00324                 }
00325             }// for yindex
00326           // end kt loop
00327           break ;
00328         } // if kt  
00329     } // for kt
00330 }
00331 //____________________________
00332 void YKPLongitudinal::AddMixedPair(const StHbtPair* pair)
00333 {
00334   // get the pair momenta 
00335   double mKt = pair->kT() ; 
00336   double mRap = pair->rap() ;
00337   double mQinv = fabs(pair->qInv()); 
00338   // caculate coulomb weight
00339   double weight=1.0;
00340   if (mCorrection)
00341     {
00342       weight = mCorrection->CoulombCorrect(pair);
00343     }
00344   
00345   // sort it into the approriate histo
00346   for(int ktindex = 0 ; ktindex < mNumberKtBins ; ktindex++)
00347     { 
00348       if( mKt >= mktBinsMin[ktindex] && mKt < mktBinsMax[ktindex] )
00349         {
00350           for(int yindex = 0 ; yindex < mNumberYBins ; yindex++)
00351             { 
00352               if( mRap >= mYBinsMin[yindex] && mRap < mYBinsMax[yindex] ) 
00353                 {
00354                   // sum up qinv pairs to do later the normalisation 
00355                   if ((mQinv < mQinvNormHi) && (mQinv > mQinvNormLo)) mNumMixedNorm[ktindex+yindex*mNumberKtBins]++;
00356                   // fill 1d qinv histo
00357                   mQinvDenominator[ktindex+yindex*mNumberKtBins].Fill(mQinv,weight) ;
00358 
00359                   // fill 3d : choose correlationfunction type and frame
00360                   if (mCtype == "YKP")
00361                     {
00362                       double qlong, qperp, q0 ;
00363                       if(mFrame =="LCMS") { pair->qYKPLCMS(qlong,qperp,q0) ;}
00364                       else if(mFrame =="CMS") { pair->qYKPCMS(qlong,qperp,q0) ;}
00365                       else if(mFrame =="PRF") { pair->qYKPPF(qlong,qperp,q0) ;}
00366                       // fill 3d histo
00367                       mDenominator[ktindex+yindex*mNumberKtBins].Fill(qlong,qperp,q0,weight) ;
00368                     }
00369                   else if (mCtype == "BP")
00370                     {
00371                       double qside, qout, qlong;
00372                       if(mFrame =="LCMS")     { qside = pair->qSideCMS(); qout = pair->qOutCMS(); qlong = pair->qLongCMS(); }
00373                       else if(mFrame =="CMS") { qside = pair->qSideBf();  qout = pair->qOutBf();  qlong = pair->qLongBf();  }
00374                       else if(mFrame =="PRF") { qside = pair->qSidePf();  qout = pair->qOutPf();  qlong = pair->qLongPf(); }
00375                       // fill 3d histo
00376                       mDenominator[ktindex+yindex*mNumberKtBins].Fill(fabs(qside),fabs(qout),fabs(qlong),weight) ;
00377                     }
00378                   // end yindex loop
00379                   break ; 
00380                 }
00381             }// for yindex
00382           // end kt loop
00383           break ;
00384         } // if kt  
00385     } // for kt
00386 }