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