AliFemtoBPLCMS3DCorrFctn.cxx

   1 ///////////////////////////////////////////////////////////////////////////
   2 //                                                                       //
   3 // AliFemtoBPLCMS3DCorrFctn: a class to calculate 3D correlation         //
   4 // for pairs of identical particles.                                     //
   5 // It also stored the weighted qinv per bin histogram for the coulomb    //
   6 // correction.                                                           //
   7 // In analysis the function should be first created in a macro, then     //
   8 // added to the analysis, and at the end of the macro the procedure to   //
   9 // write out histograms should be called.                                //
  10 //                                                                       //
  11 ///////////////////////////////////////////////////////////////////////////
  12
  13 #include "AliFemtoBPLCMS3DCorrFctn.h"
  14 //#include "AliFemtoHisto.h"
  15 #include <cstdio>
  16
  17 #ifdef __ROOT__
  18 ClassImp(AliFemtoBPLCMS3DCorrFctn)
  19 #endif
  20
  21 //____________________________
  22 AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi)
  23   :
  24   fIDNumHisto(0),
  25   fIDDenHisto(0),
  26   fIDRatHisto(0),
  27   fSMNumHisto(0),
  28   fSMDenHisto(0),
  29   fSMRatHisto(0),
  30   fCorrectionHisto(0),
  31   fCorrCFHisto(0),
  32   fNumerator(0),
  33   fDenominator(0),
  34   fRatio(0),
  35   fQinvHisto(0),
  36   fLambda(0),
  37   fRout2(0),
  38   fRside2(0),
  39   fRlong2(0),
  40   fPairCut(0),
  41   fQinvNormLo(0),
  42   fQinvNormHi(0),
  43   fNumRealsNorm(0),
  44   fNumMixedNorm(0)
  45 {
  46   // Basic constructor
  47   // set some stuff...
  48   fQinvNormLo = 0.15;
  49   fQinvNormHi = 0.18;
  50   fNumRealsNorm = 0;
  51   fNumMixedNorm = 0;
  52   //  fCorrection = 0;  // pointer to Coulomb Correction object
  53
  54   fPairCut = 0; // added Sept2000 - CorrFctn-specific PairCut
  55
  56   //  fSmearPair = 0; // no resolution correction unless user sets SmearPair
  57
  58   // set up numerator
  59   char tTitNum[100] = "Num";
  60   strcat(tTitNum,title);
  61   fNumerator = new TH3D(tTitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  62   // set up denominator
  63   char tTitDen[100] = "Den";
  64   strcat(tTitDen,title);
  65   fDenominator = new TH3D(tTitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  66   // set up uncorrected denominator
  67   char tTitDenUncoul[100] = "DenNoCoul";
  68   strcat(tTitDenUncoul,title);
  69   //  fUncorrectedDenominator = new TH3D(tTitDenUncoul,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  70   // set up ratio
  71   char tTitRat[100] = "Rat";
  72   strcat(tTitRat,title);
  73   fRatio = new TH3D(tTitRat,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  74   // set up ave qInv
  75   char tTitQinv[100] = "Qinv";
  76   strcat(tTitQinv,title);
  77   fQinvHisto = new TH3D(tTitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  78
  79   // to enable error bar calculation...
  80   fNumerator->Sumw2();
  81   fDenominator->Sumw2();
  82   //  fUncorrectedDenominator->Sumw2();
  83   fRatio->Sumw2();
  84
  85   // Following histos are for the momentum resolution correction
  86   // they are filled only if a AliFemtoSmear object is plugged in
  87   // here comes the "idea" numerator and denominator and ratio...
  88   char tTitNumID[100] = "IDNum";
  89   strcat(tTitNumID,title);
  90   fIDNumHisto = new TH3D(tTitNumID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  91   char tTitDenID[100] = "IDDen";
  92   strcat(tTitDenID,title);
  93   fIDDenHisto = new TH3D(tTitDenID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  94   char tTitRatID[100] = "IDRat";
  95   strcat(tTitRatID,title);
  96   fIDRatHisto = new TH3D(tTitRatID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  97
  98   fIDNumHisto->Sumw2();
  99   fIDDenHisto->Sumw2();
 100   fIDRatHisto->Sumw2();
 101
 102   //
 103   // here comes the "smeared" numerator and denominator...
 104   char tTitNumSM[100] = "SMNum";
 105   strcat(tTitNumSM,title);
 106   fSMNumHisto = new TH3D(tTitNumSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 107   char tTitDenSM[100] = "SMDen";
 108   strcat(tTitDenSM,title);
 109   fSMDenHisto = new TH3D(tTitDenSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 110   char tTitRatSM[100] = "SMRat";
 111   strcat(tTitRatSM,title);
 112   fSMRatHisto = new TH3D(tTitRatSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 113   //
 114   fSMNumHisto->Sumw2();
 115   fSMDenHisto->Sumw2();
 116   fSMRatHisto->Sumw2();
 117   //
 118   // here comes the correction factor (which is just ratio of ideal ratio to smeared ratio)
 119   char tTitCorrection[100] = "CorrectionFactor";
 120   strcat(tTitCorrection,title);
 121   fCorrectionHisto = new TH3D(tTitCorrection,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 122   fCorrectionHisto->Sumw2();
 123   // here comes the fully corrected correlation function
 124   char tTitCorrCF[100] = "CorrectedCF";
 125   strcat(tTitCorrCF,title);
 126   fCorrCFHisto = new TH3D(tTitCorrCF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 127   fCorrCFHisto->Sumw2();
 128
 129   // user can (and should) override these defaults...
 130   fLambda = 0.6;
 131   fRout2 = 6.0*6.0;
 132   fRside2 = 6.0*6.0;
 133   fRlong2 = 7.0*7.0;
 134
 135 }
 136
 137 AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn) :
 138   fIDNumHisto(0),
 139   fIDDenHisto(0),
 140   fIDRatHisto(0),
 141   fSMNumHisto(0),
 142   fSMDenHisto(0),
 143   fSMRatHisto(0),
 144   fCorrectionHisto(0),
 145   fCorrCFHisto(0),
 146   fNumerator(0),
 147   fDenominator(0),
 148   fRatio(0),
 149   fQinvHisto(0),
 150   fLambda(0),
 151   fRout2(0),
 152   fRside2(0),
 153   fRlong2(0),
 154   fPairCut(0),
 155   fQinvNormLo(0),
 156   fQinvNormHi(0),
 157   fNumRealsNorm(0),
 158   fNumMixedNorm(0)
 159 {
 160   // Copy constructor
 161   fIDNumHisto = new TH3D(*aCorrFctn.fIDNumHisto);
 162   fIDDenHisto = new TH3D(*aCorrFctn.fIDDenHisto);
 163   fIDRatHisto = new TH3D(*aCorrFctn.fIDRatHisto);
 164   fSMNumHisto = new TH3D(*aCorrFctn.fSMNumHisto);
 165   fSMDenHisto = new TH3D(*aCorrFctn.fSMDenHisto);
 166   fSMRatHisto = new TH3D(*aCorrFctn.fSMRatHisto);
 167   fCorrectionHisto = new TH3D(*aCorrFctn.fCorrectionHisto);
 168   fCorrCFHisto = new TH3D(*aCorrFctn.fCorrCFHisto);
 169   fNumerator = new TH3D(*aCorrFctn.fNumerator);
 170   fDenominator = new TH3D(*aCorrFctn.fDenominator);
 171   fRatio = new TH3D(*aCorrFctn.fRatio);
 172   fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
 173   fLambda = aCorrFctn.fLambda;
 174   fRout2 = aCorrFctn.fRout2;
 175   fRside2 = aCorrFctn.fRside2;
 176   fRlong2 = aCorrFctn.fRlong2;
 177   fPairCut = aCorrFctn.fPairCut;
 178   fQinvNormLo = aCorrFctn.fQinvNormLo;
 179   fQinvNormHi = aCorrFctn.fQinvNormHi;
 180   fNumRealsNorm = aCorrFctn.fNumRealsNorm;
 181   fNumMixedNorm = aCorrFctn.fNumMixedNorm;
 182 }
 183 //____________________________
 184 AliFemtoBPLCMS3DCorrFctn::~AliFemtoBPLCMS3DCorrFctn(){
 185   // Destructor
 186   delete fNumerator;
 187   delete fDenominator;
 188   delete fRatio;
 189   delete fQinvHisto;
 190   delete fIDNumHisto;
 191   delete fIDDenHisto;
 192   delete fIDRatHisto;
 193   delete fSMNumHisto;
 194   delete fSMDenHisto;
 195   delete fSMRatHisto;
 196   delete fCorrectionHisto;
 197   delete fCorrCFHisto;
 198 }
 199 //_________________________
 200 AliFemtoBPLCMS3DCorrFctn& AliFemtoBPLCMS3DCorrFctn::operator=(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn)
 201 {
 202   // assignment operator
 203   if (this == &aCorrFctn)
 204     return *this;
 205   if (fIDNumHisto) delete fIDNumHisto;
 206   fIDNumHisto = new TH3D(*aCorrFctn.fIDNumHisto);
 207   if (fIDDenHisto) delete fIDDenHisto;
 208   fIDDenHisto = new TH3D(*aCorrFctn.fIDDenHisto);
 209   if (fIDRatHisto) delete fIDRatHisto;
 210   fIDRatHisto = new TH3D(*aCorrFctn.fIDRatHisto);
 211   if (fSMNumHisto) delete fSMNumHisto;
 212   fSMNumHisto = new TH3D(*aCorrFctn.fSMNumHisto);
 213   if (fSMDenHisto) delete fSMDenHisto;
 214   fSMDenHisto = new TH3D(*aCorrFctn.fSMDenHisto);
 215   if (fSMRatHisto) delete fSMRatHisto;
 216   fSMRatHisto = new TH3D(*aCorrFctn.fSMRatHisto);
 217
 218   if (fCorrectionHisto) delete fCorrectionHisto;
 219   fCorrectionHisto = new TH3D(*aCorrFctn.fCorrectionHisto);
 220   if (fCorrCFHisto) delete fCorrCFHisto;
 221   fCorrCFHisto = new TH3D(*aCorrFctn.fCorrCFHisto);
 222   if (fNumerator) delete fNumerator;
 223   fNumerator = new TH3D(*aCorrFctn.fNumerator);
 224   if (fDenominator) delete fDenominator;
 225   fDenominator = new TH3D(*aCorrFctn.fDenominator);
 226   if (fRatio) delete fRatio;
 227   fRatio = new TH3D(*aCorrFctn.fRatio);
 228   if (fQinvHisto) delete fQinvHisto;
 229   fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
 230
 231   fLambda = aCorrFctn.fLambda;
 232   fRout2 = aCorrFctn.fRout2;
 233   fRside2 = aCorrFctn.fRside2;
 234   fRlong2 = aCorrFctn.fRlong2;
 235   fPairCut = aCorrFctn.fPairCut;
 236   fQinvNormLo = aCorrFctn.fQinvNormLo;
 237   fQinvNormHi = aCorrFctn.fQinvNormHi;
 238   fNumRealsNorm = aCorrFctn.fNumRealsNorm;
 239   fNumMixedNorm = aCorrFctn.fNumMixedNorm;
 240
 241   return *this;
 242 }
 243
 244 //_________________________
 245 void AliFemtoBPLCMS3DCorrFctn::WriteOutHistos(){
 246   // Write out all histograms to file
 247   fNumerator->Write();
 248   fDenominator->Write();
 249   //  fUncorrectedDenominator->Write();
 250   fRatio->Write();
 251   fQinvHisto->Write();
 252
 253   /*
 254     if (fSmearPair){
 255     fIDNumHisto->Write();
 256     fIDDenHisto->Write();
 257     fIDRatHisto->Write();
 258     //
 259     fSMNumHisto->Write();
 260     fSMDenHisto->Write();
 261     fSMRatHisto->Write();
 262     //
 263     fCorrectionHisto->Write();
 264     fCorrCFHisto->Write();
 265     }
 266   */
 267 }
 268
 269 //_________________________
 270 void AliFemtoBPLCMS3DCorrFctn::Finish(){
 271   // here is where we should normalize, fit, etc...
 272   double tNumFact,tDenFact;
 273   if ((fNumRealsNorm !=0) && (fNumMixedNorm !=0)){
 274     tNumFact = double(fNumRealsNorm);
 275     tDenFact = double(fNumMixedNorm);
 276   }
 277   // can happen that the fNumRealsNorm and fNumMixedNorm = 0 if you do non-standard
 278   //   things like making a new CorrFctn and just setting the Numerator and Denominator
 279   //   from OTHER CorrFctns which you read in (like when doing parallel processing)
 280   else{
 281     cout << "Warning! - no normalization constants defined - I do the best I can..." << endl;
 282     int nbins = fNumerator->GetNbinsX();
 283     int half_way = nbins/2;
 284     tNumFact = fNumerator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
 285     tDenFact = fDenominator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
 286   }
 287
 288   fRatio->Divide(fNumerator,fDenominator,tDenFact,tNumFact);
 289   //  fQinvHisto->Divide(fUncorrectedDenominator);
 290   fQinvHisto->Divide(fDenominator);
 291
 292   /*
 293   // now do all the resolution correction stuff..
 294   if (fSmearPair){  // but only do it if we have been working with a SmearPair
 295   fIDRatHisto->Divide(fIDNumHisto,fIDDenHisto);
 296   fSMRatHisto->Divide(fSMNumHisto,fSMDenHisto);
 297   fCorrectionHisto->Divide(fIDRatHisto,fSMRatHisto);
 298   fCorrCFHisto->Multiply(fRatio,fCorrectionHisto);
 299   }
 300   */
 301
 302 }
 303
 304 //____________________________
 305 AliFemtoString AliFemtoBPLCMS3DCorrFctn::Report(){
 306   // Construct the report
 307   string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n";
 308   char ctemp[100];
 309   sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
 310   stemp += ctemp;
 311   sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
 312   stemp += ctemp;
 313   sprintf(ctemp,"Number of entries in ratio:\t%E\n",fRatio->GetEntries());
 314   stemp += ctemp;
 315   sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",fQinvNormLo,fQinvNormHi);
 316   stemp += ctemp;
 317   sprintf(ctemp,"Number of pairs in Normalization region was:\n");
 318   stemp += ctemp;
 319   sprintf(ctemp,"In numerator:\t%lu\t In denominator:\t%lu\n",fNumRealsNorm,fNumMixedNorm);
 320   stemp += ctemp;
 321   /*  if (fCorrection)
 322       {
 323       float radius = fCorrection->GetRadius();
 324       sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius);
 325       }
 326       else
 327       {
 328       sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n");
 329       }
 330       stemp += ctemp;
 331   */
 332
 333   if (fPairCut){
 334     sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
 335     stemp += ctemp;
 336     stemp += fPairCut->Report();
 337   }
 338   else{
 339     sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
 340     stemp += ctemp;
 341   }
 342
 343   //
 344   AliFemtoString returnThis = stemp;
 345   return returnThis;
 346 }
 347 //____________________________
 348 void AliFemtoBPLCMS3DCorrFctn::AddRealPair( AliFemtoPair* pair){
 349   // perform operations on real pairs
 350   if (fPairCut){
 351     if (!(fPairCut->Pass(pair))) return;
 352   }
 353
 354   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 355   if ((tQinv < fQinvNormHi) && (tQinv > fQinvNormLo)) fNumRealsNorm++;
 356   double qOut = fabs(pair->QOutCMS());
 357   double qSide = fabs(pair->QSideCMS());
 358   double qLong = fabs(pair->QLongCMS());
 359
 360   fNumerator->Fill(qOut,qSide,qLong);
 361 }
 362 //____________________________
 363 void AliFemtoBPLCMS3DCorrFctn::AddMixedPair( AliFemtoPair* pair){
 364   // perform operations on mixed pairs
 365   if (fPairCut){
 366     if (!(fPairCut->Pass(pair))) return;
 367   }
 368
 369   //  double CoulombWeight = (fCorrection ? fCorrection->CoulombCorrect(pair) : 1.0);
 370   double tCoulombWeight = 1.0;
 371
 372   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 373   if ((tQinv < fQinvNormHi) && (tQinv > fQinvNormLo)) fNumMixedNorm++;
 374   double qOut = fabs(pair->QOutCMS());
 375   double qSide = fabs(pair->QSideCMS());
 376   double qLong = fabs(pair->QLongCMS());
 377
 378   fDenominator->Fill(qOut,qSide,qLong,tCoulombWeight);
 379   //  fUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0);
 380   fQinvHisto->Fill(qOut,qSide,qLong,tQinv);
 381
 382   /*
 383   // now for the momentum resolution stuff...
 384   if (fSmearPair){
 385       double CorrWeight =  1.0 +
 386       fLambda*exp((-qOut*qOut*fRout2 -qSide*qSide*fRside2 -qLong*qLong*fRlong2)/0.038936366329);
 387     CorrWeight *= CoulombWeight;  // impt.
 388
 389     fIDNumHisto->Fill(qOut,qSide,qLong,CorrWeight);
 390     fIDDenHisto->Fill(qOut,qSide,qLong,CoulombWeight);
 391
 392     fSmearPair->SetUnsmearedPair(pair);
 393     double qOut_prime = fabs(fSmearPair->SmearedPair().qOutCMS());
 394     double qSide_prime = fabs(fSmearPair->SmearedPair().qSideCMS());
 395     double qLong_prime = fabs(fSmearPair->SmearedPair().qLongCMS());
 396
 397     fSMNumHisto->Fill(qOut_prime,qSide_prime,qLong_prime,CorrWeight);
 398
 399     double SmearedCoulombWeight = ( fCorrection ?
 400                                     fCorrection->CoulombCorrect(&(fSmearPair->SmearedPair())) :
 401                                     1.0);
 402
 403     fSMDenHisto->Fill(qOut_prime,qSide_prime,qLong_prime,SmearedCoulombWeight);
 404   }
 405   */
 406 }
 407
 408