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