PWG2/FEMTOSCOPY/AliFemto/CorrFctn/AliFemtoBPLCMS3DCorrFctn.cxx

   1 /***************************************************************************
   2  *
   3  * $Id$
   4  *
   5  * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
   6  ***************************************************************************
   7  *
   8  * Description: part of STAR HBT Framework: AliFemtoMaker package
   9  *   3D Bertsch-Pratt decomposition in the LCMS frame
  10  *
  11  ***************************************************************************
  12  *
  13  * $Log$
  14  * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
  15  * First version on CVS
  16  *
  17  * Revision 1.7  2003/01/31 19:20:54  magestro
  18  * Cleared up simple compiler warnings on i386_linux24
  19  *
  20  * Revision 1.6  2002/06/07 22:51:38  lisa
  21  * Widely used AliFemtoBPLCMS3DCorrFctn class now accumulates UNcorrected denominator and has a WriteOutHistos method
  22  *
  23  * Revision 1.5  2001/05/23 00:19:04  lisa
  24  * Add in Smearing classes and methods needed for momentum resolution studies and correction
  25  *
  26  * Revision 1.4  2000/10/26 19:48:49  rcwells
  27  * Added functionality for Coulomb correction of <qInv> in 3D correltions
  28  *
  29  * Revision 1.3  2000/09/14 18:36:53  lisa
  30  * Added Qinv and ExitSep pair cuts and AliFemtoBPLCMS3DCorrFctn_SIM CorrFctn
  31  *
  32  * Revision 1.2  2000/08/23 19:43:43  lisa
  33  * added alternate normalization algorithm to 3d CorrFctns in case normal one fails
  34  *
  35  * Revision 1.1  2000/08/17 20:48:39  lisa
  36  * Adding correlationfunction in LCMS frame
  37  *
  38  *
  39  **************************************************************************/
  40
  41 #include "CorrFctn/AliFemtoBPLCMS3DCorrFctn.h"
  42 //#include "Infrastructure/AliFemtoHisto.h"
  43 #include <cstdio>
  44
  45 #ifdef __ROOT__
  46 ClassImp(AliFemtoBPLCMS3DCorrFctn)
  47 #endif
  48
  49 //____________________________
  50 AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi){
  51
  52   // set some stuff...
  53   fQinvNormLo = 0.15;
  54   fQinvNormHi = 0.18;
  55   fNumRealsNorm = 0;
  56   fNumMixedNorm = 0;
  57   //  fCorrection = 0;  // pointer to Coulomb Correction object
  58
  59   fPairCut = 0; // added Sept2000 - CorrFctn-specific PairCut
  60
  61   //  fSmearPair = 0; // no resolution correction unless user sets SmearPair
  62
  63   // set up numerator
  64   char TitNum[100] = "Num";
  65   strcat(TitNum,title);
  66   fNumerator = new TH3D(TitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  67   // set up denominator
  68   char TitDen[100] = "Den";
  69   strcat(TitDen,title);
  70   fDenominator = new TH3D(TitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  71   // set up uncorrected denominator
  72   char TitDenUncoul[100] = "DenNoCoul";
  73   strcat(TitDenUncoul,title);
  74   //  fUncorrectedDenominator = new TH3D(TitDenUncoul,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  75   // set up ratio
  76   char TitRat[100] = "Rat";
  77   strcat(TitRat,title);
  78   fRatio = new TH3D(TitRat,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  79   // set up ave qInv
  80   char TitQinv[100] = "Qinv";
  81   strcat(TitQinv,title);
  82   fQinvHisto = new TH3D(TitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  83
  84   // to enable error bar calculation...
  85   fNumerator->Sumw2();
  86   fDenominator->Sumw2();
  87   //  fUncorrectedDenominator->Sumw2();
  88   fRatio->Sumw2();
  89
  90   // Following histos are for the momentum resolution correction
  91   // they are filled only if a AliFemtoSmear object is plugged in
  92   // here comes the "idea" numerator and denominator and ratio...
  93   char TitNumID[100] = "IDNum";
  94   strcat(TitNumID,title);
  95   fIDNumHisto = new TH3D(TitNumID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  96   char TitDenID[100] = "IDDen";
  97   strcat(TitDenID,title);
  98   fIDDenHisto = new TH3D(TitDenID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  99   char TitRatID[100] = "IDRat";
 100   strcat(TitRatID,title);
 101   fIDRatHisto = new TH3D(TitRatID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 102
 103   fIDNumHisto->Sumw2();
 104   fIDDenHisto->Sumw2();
 105   fIDRatHisto->Sumw2();
 106
 107   //
 108   // here comes the "smeared" numerator and denominator...
 109   char TitNumSM[100] = "SMNum";
 110   strcat(TitNumSM,title);
 111   fSMNumHisto = new TH3D(TitNumSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 112   char TitDenSM[100] = "SMDen";
 113   strcat(TitDenSM,title);
 114   fSMDenHisto = new TH3D(TitDenSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 115   char TitRatSM[100] = "SMRat";
 116   strcat(TitRatSM,title);
 117   fSMRatHisto = new TH3D(TitRatSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 118   //
 119   fSMNumHisto->Sumw2();
 120   fSMDenHisto->Sumw2();
 121   fSMRatHisto->Sumw2();
 122   //
 123   // here comes the correction factor (which is just ratio of ideal ratio to smeared ratio)
 124   char TitCorrection[100] = "CorrectionFactor";
 125   strcat(TitCorrection,title);
 126   fCorrectionHisto = new TH3D(TitCorrection,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 127   fCorrectionHisto->Sumw2();
 128   // here comes the fully corrected correlation function
 129   char TitCorrCF[100] = "CorrectedCF";
 130   strcat(TitCorrCF,title);
 131   fCorrCFHisto = new TH3D(TitCorrCF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 132   fCorrCFHisto->Sumw2();
 133
 134   // user can (and should) override these defaults...
 135   fLambda = 0.6;
 136   fRout2 = 6.0*6.0;
 137   fRside2 = 6.0*6.0;
 138   fRlong2 = 7.0*7.0;
 139
 140 }
 141
 142 //____________________________
 143 AliFemtoBPLCMS3DCorrFctn::~AliFemtoBPLCMS3DCorrFctn(){
 144   delete fNumerator;
 145   delete fDenominator;
 146   delete fRatio;
 147   delete fQinvHisto;
 148   delete fIDNumHisto;
 149   delete fIDDenHisto;
 150   delete fIDRatHisto;
 151   delete fSMNumHisto;
 152   delete fSMDenHisto;
 153   delete fSMRatHisto;
 154   delete fCorrectionHisto;
 155   delete fCorrCFHisto;
 156 }
 157
 158 //_________________________
 159 void AliFemtoBPLCMS3DCorrFctn::WriteOutHistos(){
 160
 161   fNumerator->Write();
 162   fDenominator->Write();
 163   //  fUncorrectedDenominator->Write();
 164   fRatio->Write();
 165   fQinvHisto->Write();
 166
 167   /*
 168     if (fSmearPair){
 169     fIDNumHisto->Write();
 170     fIDDenHisto->Write();
 171     fIDRatHisto->Write();
 172     //
 173     fSMNumHisto->Write();
 174     fSMDenHisto->Write();
 175     fSMRatHisto->Write();
 176     //
 177     fCorrectionHisto->Write();
 178     fCorrCFHisto->Write();
 179     }
 180   */
 181 }
 182
 183 //_________________________
 184 void AliFemtoBPLCMS3DCorrFctn::Finish(){
 185   // here is where we should normalize, fit, etc...
 186   double NumFact,DenFact;
 187   if ((fNumRealsNorm !=0) && (fNumMixedNorm !=0)){
 188     NumFact = double(fNumRealsNorm);
 189     DenFact = double(fNumMixedNorm);
 190   }
 191   // can happen that the fNumRealsNorm and fNumMixedNorm = 0 if you do non-standard
 192   //   things like making a new CorrFctn and just setting the Numerator and Denominator
 193   //   from OTHER CorrFctns which you read in (like when doing parallel processing)
 194   else{
 195     cout << "Warning! - no normalization constants defined - I do the best I can..." << endl;
 196     int nbins = fNumerator->GetNbinsX();
 197     int half_way = nbins/2;
 198     NumFact = fNumerator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
 199     DenFact = fDenominator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
 200   }
 201
 202   fRatio->Divide(fNumerator,fDenominator,DenFact,NumFact);
 203   //  fQinvHisto->Divide(fUncorrectedDenominator);
 204   fQinvHisto->Divide(fDenominator);
 205
 206   /*
 207   // now do all the resolution correction stuff..
 208   if (fSmearPair){  // but only do it if we have been working with a SmearPair
 209   fIDRatHisto->Divide(fIDNumHisto,fIDDenHisto);
 210   fSMRatHisto->Divide(fSMNumHisto,fSMDenHisto);
 211   fCorrectionHisto->Divide(fIDRatHisto,fSMRatHisto);
 212   fCorrCFHisto->Multiply(fRatio,fCorrectionHisto);
 213   }
 214   */
 215
 216 }
 217
 218 //____________________________
 219 AliFemtoString AliFemtoBPLCMS3DCorrFctn::Report(){
 220   string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n";
 221   char ctemp[100];
 222   sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
 223   stemp += ctemp;
 224   sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
 225   stemp += ctemp;
 226   sprintf(ctemp,"Number of entries in ratio:\t%E\n",fRatio->GetEntries());
 227   stemp += ctemp;
 228   sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",fQinvNormLo,fQinvNormHi);
 229   stemp += ctemp;
 230   sprintf(ctemp,"Number of pairs in Normalization region was:\n");
 231   stemp += ctemp;
 232   sprintf(ctemp,"In numerator:\t%lu\t In denominator:\t%lu\n",fNumRealsNorm,fNumMixedNorm);
 233   stemp += ctemp;
 234   /*  if (fCorrection)
 235       {
 236       float radius = fCorrection->GetRadius();
 237       sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius);
 238       }
 239       else
 240       {
 241       sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n");
 242       }
 243       stemp += ctemp;
 244   */
 245
 246   if (fPairCut){
 247     sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
 248     stemp += ctemp;
 249     stemp += fPairCut->Report();
 250   }
 251   else{
 252     sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
 253     stemp += ctemp;
 254   }
 255
 256   //
 257   AliFemtoString returnThis = stemp;
 258   return returnThis;
 259 }
 260 //____________________________
 261 void AliFemtoBPLCMS3DCorrFctn::AddRealPair(const AliFemtoPair* pair){
 262
 263   if (fPairCut){
 264     if (!(fPairCut->Pass(pair))) return;
 265   }
 266
 267   double Qinv = fabs(pair->qInv());   // note - qInv() will be negative for identical pairs...
 268   if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumRealsNorm++;
 269   double qOut = fabs(pair->qOutCMS());
 270   double qSide = fabs(pair->qSideCMS());
 271   double qLong = fabs(pair->qLongCMS());
 272
 273   fNumerator->Fill(qOut,qSide,qLong);
 274 }
 275 //____________________________
 276 void AliFemtoBPLCMS3DCorrFctn::AddMixedPair(const AliFemtoPair* pair){
 277
 278   if (fPairCut){
 279     if (!(fPairCut->Pass(pair))) return;
 280   }
 281
 282   //  double CoulombWeight = (fCorrection ? fCorrection->CoulombCorrect(pair) : 1.0);
 283   double CoulombWeight = 1.0;
 284
 285   double Qinv = fabs(pair->qInv());   // note - qInv() will be negative for identical pairs...
 286   if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumMixedNorm++;
 287   double qOut = fabs(pair->qOutCMS());
 288   double qSide = fabs(pair->qSideCMS());
 289   double qLong = fabs(pair->qLongCMS());
 290
 291   fDenominator->Fill(qOut,qSide,qLong,CoulombWeight);
 292   //  fUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0);
 293   fQinvHisto->Fill(qOut,qSide,qLong,Qinv);
 294
 295   /*
 296   // now for the momentum resolution stuff...
 297   if (fSmearPair){
 298       double CorrWeight =  1.0 +
 299       fLambda*exp((-qOut*qOut*fRout2 -qSide*qSide*fRside2 -qLong*qLong*fRlong2)/0.038936366329);
 300     CorrWeight *= CoulombWeight;  // impt.
 301
 302     fIDNumHisto->Fill(qOut,qSide,qLong,CorrWeight);
 303     fIDDenHisto->Fill(qOut,qSide,qLong,CoulombWeight);
 304
 305     fSmearPair->SetUnsmearedPair(pair);
 306     double qOut_prime = fabs(fSmearPair->SmearedPair().qOutCMS());
 307     double qSide_prime = fabs(fSmearPair->SmearedPair().qSideCMS());
 308     double qLong_prime = fabs(fSmearPair->SmearedPair().qLongCMS());
 309
 310     fSMNumHisto->Fill(qOut_prime,qSide_prime,qLong_prime,CorrWeight);
 311
 312     double SmearedCoulombWeight = ( fCorrection ?
 313                                     fCorrection->CoulombCorrect(&(fSmearPair->SmearedPair())) :
 314                                     1.0);
 315
 316     fSMDenHisto->Fill(qOut_prime,qSide_prime,qLong_prime,SmearedCoulombWeight);
 317   }
 318   */
 319 }
 320
 321