StRoot: StHbtMaker/CorrFctn/NonIdPurityCorrFctn.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: NonIdPurityCorrFctn.cxx,v 1.2 2003/01/31 19:21:09 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  *   a simple Q-invariant correlation function           
00010  *
00011  ***************************************************************************
00012  *
00013  * $Log: NonIdPurityCorrFctn.cxx,v $
00014  * Revision 1.2  2003/01/31 19:21:09  magestro
00015  * Cleared up simple compiler warnings on i386_linux24
00016  *
00017  * Revision 1.1  2002/12/12 17:02:49  kisiel
00018  * Use KStar instead of 2*KStar for non-identical particles
00019  *
00020  * Revision 1.4  2000/01/25 17:34:45  laue
00021  * I. In order to run the stand alone version of the StHbtMaker the following
00022  * changes have been done:
00023  * a) all ClassDefs and ClassImps have been put into #ifdef __ROOT__ statements
00024  * b) unnecessary includes of StMaker.h have been removed
00025  * c) the subdirectory StHbtMaker/doc/Make has been created including everything
00026  * needed for the stand alone version
00027  *
00028  * II. To reduce the amount of compiler warning
00029  * a) some variables have been type casted
00030  * b) some destructors have been declared as virtual
00031  *
00032  * Revision 1.3  1999/07/29 02:47:09  lisa
00033  * 1) add OpeningAngle correlation function 2) add StHbtMcEventReader 3) make histos in CorrFctns do errors correctly
00034  *
00035  * Revision 1.2  1999/07/06 22:33:20  lisa
00036  * Adjusted all to work in pro and new - dev itself is broken
00037  *
00038  * Revision 1.1.1.1  1999/06/29 16:02:57  lisa
00039  * Installation of StHbtMaker
00040  *
00041  **************************************************************************/
00042 
00043 #include "StHbtMaker/CorrFctn/NonIdPurityCorrFctn.h"
00044 #include <cstdio>
00045 
00046 #ifdef __ROOT__ 
00047 ClassImp(NonIdPurityCorrFctn)
00048 #endif
00049 
00050 //____________________________
00051 NonIdPurityCorrFctn::NonIdPurityCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi, int p1Type, int p2Type){
00052   // set up numerator
00053   char TitNumP[100] = "NumP";
00054   strcat(TitNumP,title);
00055   mNumP = new StHbt1DHisto(TitNumP,title,nbins,QinvLo,QinvHi);
00056   // set up denominator
00057   char TitDenP[100] = "DenP";
00058   strcat(TitDenP,title);
00059   mDenP = new StHbt1DHisto(TitDenP,title,nbins,QinvLo,QinvHi);
00060   // set up ratio
00061   char TitRatP[100] = "RatP";
00062   strcat(TitRatP,title);
00063   mRatP = new StHbt1DHisto(TitRatP,title,nbins,QinvLo,QinvHi);
00064   // set up numerator
00065   char TitNumN[100] = "NumN";
00066   strcat(TitNumN,title);
00067   mNumN = new StHbt1DHisto(TitNumN,title,nbins,QinvLo,QinvHi);
00068   // set up denominator
00069   char TitDenN[100] = "DenN";
00070   strcat(TitDenN,title);
00071   mDenN = new StHbt1DHisto(TitDenN,title,nbins,QinvLo,QinvHi);
00072   // set up ratio
00073   char TitRatN[100] = "RatN";
00074   strcat(TitRatN,title);
00075   mRatN = new StHbt1DHisto(TitRatN,title,nbins,QinvLo,QinvHi);
00076   // set up ratio
00077   char TitRat[100] = "Rat";
00078   strcat(TitRat,title);
00079   mRat = new StHbt1DHisto(TitRat,title,nbins,QinvLo,QinvHi);
00080   // set up 
00081   char TitRatNOverP[100] = "RatNOverP";
00082   strcat(TitRatNOverP,title);
00083   mRatNOverP = new StHbt1DHisto(TitRatNOverP,title,nbins,QinvLo,QinvHi);
00084 
00085 //   // set up numerator
00086 //   char TitNumPinvP[100] = "NumPinvP";
00087 //   strcat(TitNumPinvP,title);
00088 //   mNumPinvP = new StHbt1DHisto(TitNumPinvP,title,nbins,QinvLo,QinvHi);
00089 //   // set up denominator
00090 //   char TitDenPinvP[100] = "DenPinvP";
00091 //   strcat(TitDenPinvP,title);
00092 //   mDenPinvP = new StHbt1DHisto(TitDenPinvP,title,nbins,QinvLo,QinvHi);
00093 //   // set up numerator
00094 //   char TitNumPinvN[100] = "NumPinvN";
00095 //   strcat(TitNumPinvN,title);
00096 //   mNumPinvN = new StHbt1DHisto(TitNumPinvN,title,nbins,QinvLo,QinvHi);
00097 //   // set up denominator
00098 //   char TitDenPinvN[100] = "DenPinvN";
00099 //   strcat(TitDenPinvN,title);
00100 //   mDenPinvN = new StHbt1DHisto(TitDenPinvN,title,nbins,QinvLo,QinvHi);
00101 //   // set up ratio
00102 //   char TitRatPinv[100] = "RatPinv";
00103 //   strcat(TitRatPinv,title);
00104 //   mRatPinv = new StHbt1DHisto(TitRatPinv,title,nbins,QinvLo,QinvHi);
00105 
00106 //   // set up ratio
00107 //   char TitRatPinvNormP[100] = "RatPinvNormP";
00108 //   strcat(TitRatPinvNormP,title);
00109 //   mRatPinvNormP = new StHbt1DHisto(TitRatPinvNormP,title,nbins,QinvLo,QinvHi);
00110 //   // set up ratio
00111 //   char TitRatPinvNormN[100] = "RatPinvNormN";
00112 //   strcat(TitRatPinvNormN,title);
00113 //   mRatPinvNormN = new StHbt1DHisto(TitRatPinvNormN,title,nbins,QinvLo,QinvHi);
00114 //   // set up ratio
00115 //   char TitRatPinvNorm[100] = "RatPinvNorm";
00116 //   strcat(TitRatPinvNorm,title);
00117 //   mRatPinvNorm = new StHbt1DHisto(TitRatPinvNorm,title,nbins,QinvLo,QinvHi);
00118 //   // set up ratio
00119 //   char TitRatPinvNormNOverP[100] = "RatPinvNormNOverP";
00120 //   strcat(TitRatPinvNormNOverP,title);
00121 //   mRatPinvNormNOverP = new StHbt1DHisto(TitRatPinvNormNOverP,
00122 //                                      title,nbins,QinvLo,QinvHi);
00123   mp1Type = p1Type;
00124   mp2Type = p2Type;
00125   char TitPairPurity[100];
00126   strcpy(TitPairPurity,"PairPurityOut");
00127   strcat(TitPairPurity,title);
00128   mPairPurityOut = new TProfile(TitPairPurity,title,2*nbins,-QinvHi,QinvHi);
00129   strcpy(TitPairPurity,"PairPuritySide");
00130   strcat(TitPairPurity,title);
00131   mPairPuritySide = new TProfile(TitPairPurity,title,2*nbins,-QinvHi,QinvHi);
00132   strcpy(TitPairPurity,"PairPurityLong");
00133   strcat(TitPairPurity,title);
00134   mPairPurityLong = new TProfile(TitPairPurity,title,2*nbins,-QinvHi,QinvHi);
00135 
00136   // to enable error bar calculation...
00137   mNumP->Sumw2();
00138   mDenP->Sumw2();  
00139   mRatP->Sumw2();  
00140   mNumN->Sumw2();
00141   mDenN->Sumw2();  
00142   mRatN->Sumw2();
00143   mRat->Sumw2(); 
00144   mRatNOverP->Sumw2();
00145 
00146 //   mNumPinvP->Sumw2();
00147 //   mDenPinvP->Sumw2();  
00148 //   mNumPinvN->Sumw2();
00149 //   mDenPinvN->Sumw2();  
00150 //   mRatPinv->Sumw2(); 
00151 
00152 //   mRatPinvNormP->Sumw2();
00153 //   mRatPinvNormN->Sumw2();
00154 //   mRatPinvNorm->Sumw2();
00155 //   mRatPinvNormNOverP->Sumw2();
00156 
00157 //   mQinvPt1 = new  StHbt2DHisto("QinvPt1","QinvPt1",100,0.,1.,
00158 //                             nbins,QinvLo,QinvHi);
00159 //   mQinvY1 = new  StHbt2DHisto("QinvY1","QinvY1",100,-1.,1.,
00160 //                             nbins,QinvLo,QinvHi);
00161 
00162 //   mHKCVKSame = new  StHbt2DHisto("HKCVKSame","HKCVKSame",
00163 //                               100,-1.,1.,nbins,QinvLo,QinvHi);
00164 //   mHKCVKDiff = new  StHbt2DHisto("HKCVKDiff","HKCVKDiff",
00165 //                               100,-1.,1.,nbins,QinvLo,QinvHi);
00166 }
00167 
00168 //____________________________
00169 NonIdPurityCorrFctn::~NonIdPurityCorrFctn(){
00170   delete mNumP;
00171   delete mDenP;  
00172   delete mRatP;  
00173   delete mNumN;
00174   delete mDenN;  
00175   delete mRatN;
00176   delete mRat; 
00177   delete mRatNOverP;
00178 
00179 //   delete mNumPinvP;
00180 //   delete mDenPinvP;  
00181 //   delete mNumPinvN;
00182 //   delete mDenPinvN;  
00183 //   delete mRatPinv; 
00184 
00185 //   delete mRatPinvNormP;
00186 //   delete mRatPinvNormN;
00187 //   delete mRatPinvNorm;
00188 //   delete mRatPinvNormNOverP;
00189 
00190 //   delete mQinvPt1;
00191 //   delete mQinvY1;
00192 
00193 //   delete  mHKCVKSame;
00194 //   delete  mHKCVKDiff;
00195 
00196   delete mPairPurityOut;
00197   delete mPairPuritySide;
00198   delete mPairPurityLong;
00199   
00200 }
00201 //_________________________
00202 void NonIdPurityCorrFctn::Finish(){
00203   double tScale;
00204   int tLastNormBin = mNumP->GetNbinsX();
00205   int tFirstNormBin = tLastNormBin/2+1;
00206   // Make cvk dependant correlation function
00207   mRatP->Divide(mNumP,mDenP,1.0,1.0);
00208   tScale = mRatP->Integral(tFirstNormBin,tLastNormBin);
00209   tScale/= (tLastNormBin-tFirstNormBin+1);
00210   mRatP->Scale(1./tScale);
00211   mRatN->Divide(mNumN,mDenN,1.0,1.0);
00212   tScale = mRatN->Integral(tFirstNormBin,tLastNormBin);
00213   tScale/= (tLastNormBin-tFirstNormBin+1);
00214   mRatN->Scale(1./tScale);
00215   mRatNOverP->Divide(mRatN,mRatP,1.0,1.0);
00216 
00217 //   mRatPinvNormP->Divide(mNumP,mNumPinvP,1.0,1.0);
00218 //   mRatPinvNormN->Divide(mNumN,mNumPinvN,1.0,1.0);
00219 //   mRatPinvNormNOverP->Divide(mRatPinvNormN,mRatPinvNormP);
00220 
00221   // Regular correlation function
00222   TH1D tHNum(*mNumP);
00223   tHNum.SetName("tHNum");
00224   tHNum.Add(mNumN);
00225   TH1D tHDen(*mDenP);
00226   tHDen.SetName("tHDen");
00227   tHDen.Add(mDenN); 
00228   mRat->Divide(&tHNum,&tHDen);
00229   tScale = mRat->Integral(tFirstNormBin,tLastNormBin);
00230   tScale/= (tLastNormBin-tFirstNormBin+1);
00231   mRat->Scale(1./tScale);
00232 
00233   // Pinv correlation function
00234 //   TH1D tHNumPinv(*mNumPinvP);
00235 //   tHNumPinv.SetName("tHNumPinv");
00236 //   tHNumPinv.Add(mNumPinvN);
00237 //   TH1D tHDenPinv(*mDenPinvP);
00238 //   tHDenPinv.SetName("tHDenPinv");
00239 //   tHDenPinv.Add(mDenPinvN); 
00240 //   mRatPinv->Divide(&tHNumPinv,&tHDenPinv);
00241 //   tScale = mRatPinv->Integral(tFirstNormBin,tLastNormBin);
00242 //   tScale/= (tLastNormBin-tFirstNormBin+1);
00243 //   mRatPinv->Scale(1./tScale);
00244 
00245   // Pinv normalised correlation function
00246 //   mRatPinvNorm->Divide(&tHNum,&tHNumPinv);
00247 }
00248 
00249 void NonIdPurityCorrFctn::Write(){
00250   mNumP->Write();
00251   mDenP->Write();  
00252   mRatP->Write();  
00253   mNumN->Write();
00254   mDenN->Write();  
00255   mRatN->Write();
00256   mRat->Write(); 
00257   mRatNOverP->Write();
00258 
00259 //   mNumPinvP->Write();
00260 //   mDenPinvP->Write();  
00261 //   mNumPinvN->Write();
00262 //   mDenPinvN->Write();  
00263 //   mRatPinv->Write(); 
00264 
00265 //   mRatPinvNormP->Write();
00266 //   mRatPinvNormN->Write();
00267 //   mRatPinvNorm->Write();
00268 //   mRatPinvNormNOverP->Write();
00269 
00270 //   mHKCVKSame->Write();
00271 //   mHKCVKDiff->Write();
00272   //mQinvPt1->Write();
00273   //mQinvY1->Write();
00274 
00275   mPairPurityOut->Write();  
00276   mPairPuritySide->Write();  
00277   mPairPurityLong->Write();  
00278 }
00279 
00280 //____________________________
00281 StHbtString NonIdPurityCorrFctn::Report(){
00282   string stemp = "Qinv Correlation Function Report:\n";
00283   char ctemp[100];
00284   sprintf(ctemp,"Number of entries in numerator:\t%E\n",mNumP->GetEntries());
00285   stemp += ctemp;
00286   sprintf(ctemp,"Number of entries in denominator:\t%E\n",mDenP->GetEntries());
00287   stemp += ctemp;
00288   sprintf(ctemp,"Number of entries in ratio:\t%E\n",mRatP->GetEntries());
00289   stemp += ctemp;
00290   StHbtString returnThis = stemp;
00291   return returnThis;
00292 }
00293 //____________________________
00294 void NonIdPurityCorrFctn::AddRealPair(const StHbtPair* pair){
00295   double tKStar = fabs(pair->KStar());
00296   //double tPStar = 2*fabs(pair->KStarFlipped());
00297   double tCVK = pair->CVK();
00298   double pPurity = 0.0;
00299   double tKOut = pair->dKOut();
00300   double tKSide = pair->dKSide();
00301   double tKLong = pair->dKLong();
00302   
00303   //mQinvPt1->Fill(pair->track1()->FourMomentum().perp(),tKStar,1);
00304   //mQinvY1->Fill(pair->track1()->FourMomentum().rapidity(),tKStar,1);
00305   
00306   //if(pair->track1()->FourMomentum().pz()<0.){
00307   //  cout << pair->track1()->FourMomentum() << endl;
00308   //}
00309 //   mHKCVKSame->Fill(tCVK,tKStar,1.);
00310 
00311   if(tCVK>0.){
00312     mNumP->Fill(tKStar);
00313   }
00314   else{
00315     mNumN->Fill(tKStar);
00316   }
00317 //   if(pair->CVKFlipped()>0.){
00318 //     mNumPinvP->Fill(tPStar);
00319 //   }
00320 //   else{
00321 //     mNumPinvN->Fill(tPStar);
00322 //   }
00323   
00324   switch (mp1Type)
00325     {
00326     case 1:
00327       pPurity =  pair->track1()->GetPionPurity();
00328       break;
00329     case 2:
00330       pPurity =  pair->track1()->GetKaonPurity();
00331       break;
00332     case 3:
00333       pPurity =  pair->track1()->GetProtonPurity();
00334       break;
00335     case 4:
00336       pPurity =  pair->track1()->Track()->PidProbElectron();
00337       break;
00338     }
00339 
00340   //  cout << "First particle purity: " << pPurity << endl;
00341   
00342   switch (mp2Type)
00343     {
00344     case 1:
00345       pPurity *=  pair->track2()->GetPionPurity();
00346       break;
00347     case 2:
00348       pPurity *=  pair->track2()->GetKaonPurity();
00349       break;
00350     case 3:
00351       pPurity *=  pair->track2()->GetProtonPurity();
00352       break;
00353     case 4:
00354       pPurity *=  pair->track2()->Track()->PidProbElectron();
00355       break;
00356     }
00357   //  cout << "Pair purity: " << pPurity << endl;
00358   
00359   mPairPurityOut->Fill(((tKOut>0)*2-1) *tKStar, pPurity);
00360   mPairPuritySide->Fill(((tKSide>0)*2-1)*tKStar, pPurity);
00361   mPairPurityLong->Fill(((tKLong>0)*2-1)*tKStar, pPurity);
00362   
00363 }
00364 //____________________________
00365 void NonIdPurityCorrFctn::AddMixedPair(const StHbtPair* pair){
00366   double tKStar = 2*fabs(pair->KStar());
00367   //double tPStar = 2*fabs(pair->KStarFlipped());
00368   double tCVK = pair->CVK();
00369 //   mHKCVKDiff->Fill(tCVK,tKStar,1.);
00370   if(tCVK>0.){
00371     mDenP->Fill(tKStar);
00372   }
00373   else{
00374     mDenN->Fill(tKStar);
00375   }
00376 //   if(pair->CVKFlipped()>0.){
00377 //     mDenPinvP->Fill(tPStar);
00378 //   }
00379 //   else{
00380 //     mDenPinvN->Fill(tPStar); 
00381 //   }
00382 }
00383 
00384