PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 ///                                                                          ///
   3 /// AliFemtoModelBPLCMSCorrFctn - the class for correlation function which   ///
   4 /// uses the model framework and weight generation and calculated the 3D     ///
   5 /// correlation function in the Bertsh-Pratt LCMS system                     ///
   6 /// Authors: Adam Kisiel, kisiel@mps.ohio-state.edu                          ///
   7 ///                                                                          ///
   8 ////////////////////////////////////////////////////////////////////////////////
   9 #include "AliFemtoModelBPLCMSCorrFctn.h"
  10 #include "AliFemtoPair.h"
  11 #include "AliFemtoModelManager.h"
  12 #include "AliFemtoKTPairCut.h"
  13 #include "AliFemtoAnalysisReactionPlane.h"
  14 #include <cstdio>
  15
  16 #ifdef __ROOT__
  17 ClassImp(AliFemtoModelBPLCMSCorrFctn)
  18 #endif
  19
  20 //____________________________
  21 AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi)
  22   :
  23   AliFemtoModelCorrFctn(title, nbins, QLo, QHi),
  24   fNumerator3DTrue(0),
  25   fNumerator3DFake(0),
  26   fDenominator3D(0),
  27   fQinvHisto(0),
  28   fPairCut(0),
  29   fUseRPSelection(0)
  30 {
  31   // set up true numerator
  32   char tTitNumT[100] = "Num3DTrue";
  33   strcat(tTitNumT,title);
  34   fNumerator3DTrue = new TH3D(tTitNumT,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  35   // set up fake numerator
  36   char tTitNumF[100] = "Num3DFake";
  37   strcat(tTitNumF,title);
  38   fNumerator3DFake = new TH3D(tTitNumF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  39   // set up denominator
  40   char tTitDen[100] = "Den3D";
  41   strcat(tTitDen,title);
  42   fDenominator3D = new TH3D(tTitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  43   // set up ave qInv
  44   char tTitQinv[100] = "Qinv";
  45   strcat(tTitQinv,title);
  46   fQinvHisto = new TH3D(tTitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  47
  48   // to enable error bar calculation...
  49   fNumerator3DTrue->Sumw2();
  50   fNumerator3DFake->Sumw2();
  51   fDenominator3D->Sumw2();
  52 }
  53
  54 AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(const AliFemtoModelBPLCMSCorrFctn& aCorrFctn) :
  55   AliFemtoModelCorrFctn(aCorrFctn),
  56   fNumerator3DTrue(0),
  57   fNumerator3DFake(0),
  58   fDenominator3D(0),
  59   fQinvHisto(0),
  60   fPairCut(0),
  61   fUseRPSelection(0)
  62 {
  63   // Copy constructor
  64   fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue);
  65   fNumerator3DFake = new TH3D(*aCorrFctn.fNumerator3DFake);
  66   fDenominator3D   = new TH3D(*aCorrFctn.fDenominator3D);
  67   fQinvHisto       = new TH3D(*aCorrFctn.fQinvHisto);
  68   fPairCut         = aCorrFctn.fPairCut->Clone();
  69 }
  70 //____________________________
  71 AliFemtoModelBPLCMSCorrFctn::~AliFemtoModelBPLCMSCorrFctn()
  72 {
  73   // destructor
  74   if (fNumeratorTrue) delete fNumeratorTrue;
  75   if (fNumeratorFake) delete fNumeratorFake;
  76   if (fDenominator) delete fDenominator;
  77   delete fNumerator3DTrue;
  78   delete fNumerator3DFake;
  79   delete fDenominator3D;
  80   delete fQinvHisto;
  81   if (fPairCut) delete fPairCut;
  82 }
  83 //_________________________
  84 AliFemtoModelBPLCMSCorrFctn& AliFemtoModelBPLCMSCorrFctn::operator=(const AliFemtoModelBPLCMSCorrFctn& aCorrFctn)
  85 {
  86   // Assignment operator
  87   if (this == &aCorrFctn)
  88     return *this;
  89   if (fNumerator3DTrue) delete fNumerator3DTrue;
  90   fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue);
  91   if (fNumerator3DFake) delete fNumerator3DFake;
  92   fNumerator3DFake = new TH3D(*aCorrFctn.fNumerator3DFake);
  93   if (fDenominator3D) delete fDenominator3D;
  94   fDenominator3D = new TH3D(*aCorrFctn.fDenominator3D);
  95   if (fQinvHisto) delete fQinvHisto;
  96   fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
  97   fPairCut = aCorrFctn.fPairCut->Clone();
  98
  99   return *this;
 100 }
 101
 102 //_________________________
 103 void AliFemtoModelBPLCMSCorrFctn::Write(){
 104   // Write out data histograms
 105   AliFemtoModelCorrFctn::Write();
 106   fNumerator3DTrue->Write();
 107   fNumerator3DFake->Write();
 108   fDenominator3D->Write();
 109   fQinvHisto->Write();
 110 }
 111 //________________________
 112 TList* AliFemtoModelBPLCMSCorrFctn::GetOutputList()
 113 {
 114   // Prepare the list of objects to be written to the output
 115   TList *tOutputList = AliFemtoModelCorrFctn::GetOutputList();
 116
 117   tOutputList->Add(fNumerator3DTrue);
 118   tOutputList->Add(fNumerator3DFake);
 119   tOutputList->Add(fDenominator3D);
 120   tOutputList->Add(fQinvHisto);
 121
 122   return tOutputList;
 123 }
 124
 125 //_________________________
 126 void AliFemtoModelBPLCMSCorrFctn::Finish(){
 127   fQinvHisto->Divide(fDenominator);
 128 }
 129
 130 //____________________________
 131 AliFemtoString AliFemtoModelBPLCMSCorrFctn::Report(){
 132   // Prepare a report from the execution
 133   string stemp = "LCMS Frame Bertsch-Pratt 3D Model Correlation Function Report:\n";
 134   char ctemp[100];
 135   sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumeratorTrue->GetEntries());
 136   stemp += ctemp;
 137   sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
 138   stemp += ctemp;
 139   /*  if (fCorrection)
 140       {
 141       float radius = fCorrection->GetRadius();
 142       sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius);
 143       }
 144       else
 145       {
 146       sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n");
 147       }
 148       stemp += ctemp;
 149   */
 150
 151   //
 152   AliFemtoString returnThis = stemp;
 153   return returnThis;
 154 }
 155 //____________________________
 156 void AliFemtoModelBPLCMSCorrFctn::AddRealPair( AliFemtoPair* pair)
 157 {
 158   // Store a real pair in numerator
 159   if (fPairCut){
 160     if (fUseRPSelection) {
 161       AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
 162       if (!ktc) {
 163         cout << "RP aware cut requested, but not connected to the CF" << endl;
 164         if (!(fPairCut->Pass(pair))) return;
 165       }
 166       else {
 167         AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
 168         if (!arp) {
 169           cout << "RP aware cut requested, but not connected to the CF" << endl;
 170           if (!(fPairCut->Pass(pair))) return;
 171         }
 172         if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
 173       }
 174     }
 175     else
 176       if (!(fPairCut->Pass(pair))) return;
 177   }
 178 //   if (fPairCut){
 179 //     if (!(fPairCut->Pass(pair))) return;
 180 //   }
 181
 182   Double_t weight = fManager->GetWeight(pair);
 183
 184   double qOut = (pair->QOutCMS());
 185   double qSide = (pair->QSideCMS());
 186   double qLong = (pair->QLongCMS());
 187
 188   fNumerator3DTrue->Fill(qOut, qSide, qLong, weight);
 189   fNumeratorTrue->Fill(pair->QInv(), weight);
 190 }
 191 //____________________________
 192 void AliFemtoModelBPLCMSCorrFctn::AddMixedPair( AliFemtoPair* pair){
 193   // store mixed pair in denominator
 194   if (fPairCut){
 195     if (fUseRPSelection) {
 196       AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
 197       if (!ktc) {
 198         cout << "RP aware cut requested, but not connected to the CF" << endl;
 199         if (!(fPairCut->Pass(pair))) return;
 200       }
 201       else {
 202         AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
 203         if (!arp) {
 204           cout << "RP aware cut requested, but not connected to the CF" << endl;
 205           if (!(fPairCut->Pass(pair))) return;
 206         }
 207         if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
 208       }
 209     }
 210     else
 211       if (!(fPairCut->Pass(pair))) return;
 212   }
 213 //   if (fPairCut){
 214 //     if (!(fPairCut->Pass(pair))) return;
 215 //   }
 216
 217   Double_t weight = fManager->GetWeight(pair);
 218
 219   double qOut = (pair->QOutCMS());
 220   double qSide = (pair->QSideCMS());
 221   double qLong = (pair->QLongCMS());
 222
 223   fNumerator3DFake->Fill(qOut, qSide, qLong, weight);
 224   fDenominator3D->Fill(qOut, qSide, qLong, 1.0);
 225   fNumeratorFake->Fill(pair->QInv(), weight);
 226   fDenominator->Fill(pair->QInv(), 1.0);
 227
 228 }
 229 //_______________________
 230 AliFemtoModelCorrFctn* AliFemtoModelBPLCMSCorrFctn::Clone()
 231 {
 232   // Clone the correlation function
 233   AliFemtoModelBPLCMSCorrFctn *tCopy = new AliFemtoModelBPLCMSCorrFctn(*this);
 234
 235   return tCopy;
 236 }
 237
 238 void AliFemtoModelBPLCMSCorrFctn::SetSpecificPairCut(AliFemtoPairCut* aCut)
 239 {
 240   fPairCut = aCut;
 241 }
 242
 243 void AliFemtoModelBPLCMSCorrFctn::SetUseRPSelection(unsigned short aRPSel)
 244 {
 245   fUseRPSelection = aRPSel;
 246 }