StRoot: StHbtMaker/CorrFctn/Q3invCorrFctn.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: Q3invCorrFctn.cxx,v 1.4 2001/06/03 21:05:42 willson Exp $
00004  *
00005  * Author: Robert Willson, Ohio State, willson@bnl.gov
00006  ***************************************************************************
00007  *  
00008  * Description: part of STAR HBT Framework: StHbtMaker package
00009  *   A simple Q-invariant correlation function for three particle analyses.    
00010  * 
00011  ***************************************************************************
00012  *
00013  * $Log: Q3invCorrFctn.cxx,v $
00014  * Revision 1.4  2001/06/03 21:05:42  willson
00015  * Bins in entrance separation
00016  *
00017  * Revision 1.2  2000/04/12 01:53:28  willson
00018  * Initial Installation - Comments Added
00019  *
00020  * 
00021  ***************************************************************************/
00022 
00023 #include "StHbtMaker/CorrFctn/Q3invCorrFctn.h"
00024 //#include "StHbtMaker/Infrastructure/StHbtHisto.hh"
00025 #include <cstdio>
00026 
00027 ClassImp(Q3invCorrFctn)
00028 
00029 //____________________________
00030 Q3invCorrFctn::Q3invCorrFctn(char* title, const int& nbinsQ, const float& QinvLo, const float& QinvHi, const int& nbinsMerge, const float& MergeLo, const float& MergeHi, const float& Split){
00031   // This one does not do a quality calulation
00032   mSplit = Split;
00033   mPHist = 0;
00034   // set up numerator
00035   //  title = "Num Q3inv (MeV/c)";
00036   char TitNum[100] = "Num";
00037   strcat(TitNum,title);
00038   mNumerator = new StHbt2DHisto(TitNum,title,nbinsQ,QinvLo,QinvHi,nbinsMerge,MergeLo,MergeHi);
00039   // set up denominator
00040   //title = "Den Q3inv (MeV/c)";
00041   char TitDen[100] = "Den";
00042   strcat(TitDen,title);
00043   mDenominator = new StHbt2DHisto(TitDen,title,nbinsQ,QinvLo,QinvHi,nbinsMerge,MergeLo,MergeHi);
00044   // set up ratio
00045   //title = "Ratio Q3inv (MeV/c)";
00046   char TitRat[100] = "Rat";
00047   strcat(TitRat,title);
00048   mRatio = new StHbt2DHisto(TitRat,title,nbinsQ,QinvLo,QinvHi,nbinsMerge,MergeLo,MergeHi);
00049   // this next bit is unfortunately needed so that we can have many histos of same "title"
00050   // it is neccessary if we typedef StHbt1DHisto to TH1d (which we do)
00051   //mNumer->SetDirectory(0);
00052   //mDenom->SetDirectory(0);
00053   //mRatio->SetDirectory(0);
00054 
00055   // to enable error bar calculation...
00056   mNumerator->Sumw2();
00057   mDenominator->Sumw2();
00058   mRatio->Sumw2();
00059 
00060 }
00061 
00062 //____________________________
00063 Q3invCorrFctn::~Q3invCorrFctn(){
00064   delete mNumerator;
00065   delete mDenominator;
00066   delete mRatio;
00067 }
00068 
00069 //_________________________
00070 void Q3invCorrFctn::Finish(){
00071   mRatio->Divide(mNumerator,mDenominator,1.0,1.0);    
00072 }
00073 
00074 //____________________________
00075 StHbtString Q3invCorrFctn::Report(){
00076   string stemp = "Q3inv Correlation Function Report:\n";
00077   char ctemp[100];
00078   sprintf(ctemp,"Number of entries in numerator:\t%E\n",mNumerator->GetEntries());
00079   stemp += ctemp;
00080   sprintf(ctemp,"Number of entries in denominator:\t%E\n",mDenominator->GetEntries());
00081   stemp += ctemp;
00082   sprintf(ctemp,"Number of entries in ratio:\t%E\n",mRatio->GetEntries());
00083   stemp += ctemp;
00084   //  stemp += mCoulombWeight->Report();
00085   StHbtString returnThis = stemp;
00086   return returnThis;
00087 }
00088 //____________________________
00089 void Q3invCorrFctn::AddRealTriplet(const StHbtTriplet* triplet){
00090   if (mPHist) {
00091     mPHist->Fill(triplet->track1()->Track()->P().mag());
00092     mPHist->Fill(triplet->track2()->Track()->P().mag());
00093     mPHist->Fill(triplet->track3()->Track()->P().mag());
00094   }
00095   double Qual = triplet->quality();
00096   double entSep = triplet->NominalTpcEntranceSeparation();
00097   double Q3inv = fabs(triplet->qInv());   // note - qInv() will be negative for identical triplets...  
00098   if (Qual<mSplit) mNumerator->Fill(Q3inv, entSep, 1.0);
00099 }
00100 
00101 //____________________________
00102 void Q3invCorrFctn::AddMixedTriplet(const StHbtTriplet* triplet){
00103   double weight = 1.0;
00104   if (mCorrection.GetRadius()!=-1) {
00105     StHbtPair pair1(triplet->track1(), triplet->track2());
00106     StHbtPair pair2(triplet->track2(), triplet->track3());
00107     StHbtPair pair3(triplet->track3(), triplet->track1());
00108     
00109     double weight1 = mCorrection.CoulombCorrect(&pair1);
00110     double weight2 = mCorrection.CoulombCorrect(&pair2);
00111     double weight3 = mCorrection.CoulombCorrect(&pair3);
00112     
00113     weight = weight1*weight2*weight3;
00114   }
00115   else 
00116     weight = 1.0;
00117   
00118   double Qual = triplet->quality();
00119   double entSep = triplet->NominalTpcEntranceSeparation();
00120   double Q3inv = fabs(triplet->qInv());   // note - qInv() will be negative for identical triplets...  
00121   if (Qual<mSplit) mDenominator->Fill(Q3inv, entSep, weight);
00122 }
00123 //____________________________
00124 void Q3invCorrFctn::AddCorrection(StHbtCoulomb* correction) {
00125   mCorrection = *correction;
00126 }