StRoot: StHbtMaker/ThCorrFctn/ThBPCorrFctn.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: 
00004  *
00005  * Author: Laurent Conin, Fabrice Retiere, Subatech, France
00006  ***************************************************************************
00007  *
00008  * Description : Calculate the theoretical QInv correlation function 
00009  *
00010  ***************************************************************************
00011  *
00012  * $Log: 
00013  *
00014  ***************************************************************************/
00015 
00016 #include "Stsstream.h"
00017 #include "StHbtMaker/ThCorrFctn/ThBPCorrFctn.h"
00018 
00019 #ifdef __ROOT__
00020 ClassImp(ThBPCorrFctn)
00021 #endif
00022 
00023 //____________________________
00024 void ThBPCorrFctn::AddNum(StHbtThPair* aThPair){
00025 
00026   double tQinv = 2 * fabs(aThPair->GetMeasPair()->KStar());   // note - qInv() will be negative for identical pairs...
00027   if ((tQinv < mQinvNormHi) && (tQinv > mQinvNormLo)) mNumRealsNorm++;
00028   double qOut = aThPair->GetMeasPair()->qOutCMS();
00029   double qSide = aThPair->GetMeasPair()->qSideCMS();
00030   double qLong = fabs(aThPair->GetMeasPair()->qLongCMS());
00031 
00032   //  cout << "Components Num: inv out side long: " << tQinv << " " << qOut << " " << qSide << " " << qLong << endl;
00033   mNumerator->Fill(qOut,qSide,qLong,aThPair->GetWeightNum());
00034 }
00035 void ThBPCorrFctn::AddDen(StHbtThPair* aThPair){
00036 
00037   double tQinv= 2 * fabs(aThPair->GetMeasPair()->KStar());
00038   if ((tQinv < mQinvNormHi) && (tQinv > mQinvNormLo)) mNumMixedNorm++;
00039   double qOut = aThPair->GetMeasPair()->qOutCMS();
00040   double qSide = aThPair->GetMeasPair()->qSideCMS();
00041   double qLong = fabs(aThPair->GetMeasPair()->qLongCMS());
00042 
00043   if ((addedHistos == 1) && (tQinv < 0.1)) {
00044     add2DHistos[0]->Fill(fabs(qOut), -(fabs(qOut) - fabs(aThPair->RealqOutCMS())), 1.0);
00045     add2DHistos[1]->Fill(fabs(qSide), -(fabs(qSide) - fabs(aThPair->RealqSideCMS())), 1.0);
00046     add2DHistos[2]->Fill(fabs(qLong), -(fabs(qLong) - fabs(aThPair->RealqLongCMS())), 1.0);
00047   }
00048 
00049   //  cout << "Components Den: inv out side long: " << tQinv << " " << qOut << " " << qSide << " " << qLong << endl;
00050   mDenominator->Fill(qOut,qSide,qLong,aThPair->GetWeightDen());
00051   mQinvHisto->Fill(qOut,qSide,qLong,tQinv);
00052 }
00053 
00054 void ThBPCorrFctn::Finish(){
00055   // here is where we should normalize, fit, etc...
00056   double NumFact,DenFact;
00057   if ((mNumRealsNorm !=0) && (mNumMixedNorm !=0)){
00058     NumFact = double(mNumRealsNorm);
00059     DenFact = double(mNumMixedNorm);
00060   }
00061   // can happen that the mNumRealsNorm and mNumMixedNorm = 0 if you do non-standard
00062   //   things like making a new CorrFctn and just setting the Numerator and Denominator
00063   //   from OTHER CorrFctns which you read in (like when doing parallel processing) 
00064   else{
00065     cout << "Warning! - no normalization constants defined - I do the best I can..." << endl;
00066     int nbins = mNumerator->GetNbinsX();
00067     int half_way = nbins/2;
00068     NumFact = mNumerator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
00069     DenFact = mDenominator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
00070   }
00071   
00072   mRatio->Divide(mNumerator,mDenominator,DenFact,NumFact);
00073   mQinvHisto->Divide(mDenominator);
00074 }
00075 
00076 ThBPCorrFctn::ThBPCorrFctn(char* aTitle, int aNBins, 
00077                            double aHLo, double aHHi, int addHistos)
00078 {
00079   // set some stuff...
00080   mQinvNormLo = 0.15;
00081   mQinvNormHi = 0.18;
00082   mNumRealsNorm = 0;
00083   mNumMixedNorm = 0;
00084 
00085   // set up numerator
00086   char TitNum[100] = "Num";
00087   strcat(TitNum,aTitle);
00088   mNumerator = new StHbt3DHisto(TitNum,aTitle,aNBins,aHLo,aHHi,aNBins,aHLo,aHHi,aNBins/2,0.0,aHHi);
00089   // set up denominator
00090   char TitDen[100] = "Den";
00091   strcat(TitDen,aTitle);
00092   mDenominator = new StHbt3DHisto(TitDen,aTitle,aNBins,aHLo,aHHi,aNBins,aHLo,aHHi,aNBins/2,0.0,aHHi);
00093   // set up ratio
00094   char TitRat[100] = "Rat";
00095   strcat(TitRat,aTitle);
00096   mRatio = new StHbt3DHisto(TitRat,aTitle,aNBins,aHLo,aHHi,aNBins,aHLo,aHHi,aNBins/2,0.0,aHHi);
00097   // set up ave qInv
00098   char TitQinv[100] = "Qinv";
00099   strcat(TitQinv,aTitle);
00100   mQinvHisto = new StHbt3DHisto(TitQinv,aTitle,aNBins,aHLo,aHHi,aNBins,aHLo,aHHi,aNBins/2,0.0,aHHi);
00101 
00102   // to enable error bar calculation...
00103   mNumerator->Sumw2();
00104   mDenominator->Sumw2();
00105   mRatio->Sumw2();
00106   
00107   addedHistos = addHistos;
00108   if (addHistos == 1) {
00109     Int_t qdepNBins = 10;
00110     Int_t qdiffNBins = 50;
00111     Float_t qdepMin = 0.0, qdepMax=0.1;
00112     Float_t qdiffMin = -0.025, qdiffMax = 0.025;
00113     numAdd2DHistos = 3;
00114     typedef StHbt2DHisto* StHbt2DHistoP;
00115     add2DHistos = new StHbt2DHistoP[3];
00116     char Tit1add[100] = "Qout";
00117     strcat(Tit1add,aTitle);
00118     add2DHistos[0] = new StHbt2DHisto(Tit1add,aTitle,qdepNBins, qdepMin, qdepMax, qdiffNBins, qdiffMin, qdiffMax);
00119     char Tit2add[100] = "Qside";
00120     strcat(Tit2add,aTitle);
00121     add2DHistos[1] = new StHbt2DHisto(Tit2add,aTitle,qdepNBins, qdepMin, qdepMax, qdiffNBins, qdiffMin, qdiffMax);
00122     char Tit3add[100] = "Qlong";
00123     strcat(Tit3add,aTitle);
00124     add2DHistos[2] = new StHbt2DHisto(Tit3add,aTitle,qdepNBins, qdepMin, qdepMax, qdiffNBins, qdiffMin, qdiffMax);
00125   }
00126 }
00127 
00128 ThBPCorrFctn::ThBPCorrFctn(const ThBPCorrFctn& ThCf)
00129 {
00130   mQinvNormLo = ThCf.mQinvNormLo;
00131   mQinvNormHi = ThCf.mQinvNormHi;
00132   mNumRealsNorm = ThCf.mNumRealsNorm;
00133   mNumMixedNorm = ThCf.mNumMixedNorm;
00134   
00135   mNumerator = new StHbt3DHisto(*ThCf.mNumerator);
00136   mDenominator = new StHbt3DHisto(*ThCf.mDenominator);
00137   mRatio = new StHbt3DHisto(*ThCf.mRatio);
00138   mQinvHisto = new StHbt3DHisto(*ThCf.mQinvHisto);
00139 }
00140 
00141 StHbtString ThBPCorrFctn::Report() 
00142 {
00143   ostrstream *out = new ostrstream();
00144   (*out) << "Report form the ThBPLCMS Correlation function" << endl;
00145   return out->str();
00146 }
00147 
00148 ThBPCorrFctn::~ThBPCorrFctn() 
00149 {
00150   delete mNumerator;
00151   delete mDenominator;
00152   delete mRatio;
00153   delete mQinvHisto;
00154   if (addedHistos > 0) {
00155     for (int i=0; i<numAdd2DHistos; i++)
00156       delete add2DHistos[i];
00157     delete add2DHistos;
00158   }
00159 }
00160 
00161 void ThBPCorrFctn::Write()  {
00162   mNumerator->Write();mDenominator->Write();mRatio->Write();
00163   if (addedHistos > 0) {
00164     for (int i=0; i<numAdd2DHistos; i++)
00165       add2DHistos[i]->Write();
00166   }
00167 };
00168 
00169 inline StHbt3DHisto* ThBPCorrFctn::Numerator() const {  return mNumerator;};
00170 inline StHbt3DHisto* ThBPCorrFctn::Denominator() const { return mDenominator;};
00171 inline StHbt3DHisto* ThBPCorrFctn::Ratio() const { return mRatio;};
00172 inline StHbtThCorrFctn* ThBPCorrFctn::ThClone() const {return new ThBPCorrFctn(*this);}