PWG2/FEMTOSCOPY/AliFemto/AliFemtoQinvCorrFctn.cxx

   1 ///////////////////////////////////////////////////////////////////////////
   2 //                                                                       //
   3 // AliFemtoQinvCorrFctn:                                                 //
   4 // a simple Q-invariant correlation function                             //
   5 //                                                                       //
   6 ///////////////////////////////////////////////////////////////////////////
   7
   8 #include "AliFemtoQinvCorrFctn.h"
   9 //#include "AliFemtoHisto.h"
  10 #include <cstdio>
  11
  12 #ifdef __ROOT__
  13 ClassImp(AliFemtoQinvCorrFctn)
  14 #endif
  15
  16 //____________________________
  17 AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
  18   fNumerator(0),
  19   fDenominator(0),
  20   fRatio(0),
  21   fkTMonitor(0)
  22 {
  23   // set up numerator
  24   //  title = "Num Qinv (MeV/c)";
  25   char tTitNum[101] = "Num";
  26   strncat(tTitNum,title, 100);
  27   fNumerator = new TH1D(tTitNum,title,nbins,QinvLo,QinvHi);
  28   // set up denominator
  29   //title = "Den Qinv (MeV/c)";
  30   char tTitDen[101] = "Den";
  31   strncat(tTitDen,title, 100);
  32   fDenominator = new TH1D(tTitDen,title,nbins,QinvLo,QinvHi);
  33   // set up ratio
  34   //title = "Ratio Qinv (MeV/c)";
  35   char tTitRat[101] = "Rat";
  36   strncat(tTitRat,title, 100);
  37   fRatio = new TH1D(tTitRat,title,nbins,QinvLo,QinvHi);
  38   char tTitkT[101] = "kTDep";
  39   strncat(tTitkT,title, 100);
  40   fkTMonitor = new TH1D(tTitkT,title,200,0.0,2.0);
  41   // this next bit is unfortunately needed so that we can have many histos of same "title"
  42   // it is neccessary if we typedef TH1D to TH1d (which we do)
  43   //fNumerator->SetDirectory(0);
  44   //fDenominator->SetDirectory(0);
  45   //fRatio->SetDirectory(0);
  46
  47   // to enable error bar calculation...
  48   fNumerator->Sumw2();
  49   fDenominator->Sumw2();
  50   fRatio->Sumw2();
  51   fkTMonitor->Sumw2();
  52 }
  53
  54 //____________________________
  55 AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(const AliFemtoQinvCorrFctn& aCorrFctn) :
  56   AliFemtoCorrFctn(),
  57   fNumerator(0),
  58   fDenominator(0),
  59   fRatio(0),
  60   fkTMonitor(0)
  61 {
  62   // copy constructor
  63   fNumerator = new TH1D(*aCorrFctn.fNumerator);
  64   fDenominator = new TH1D(*aCorrFctn.fDenominator);
  65   fRatio = new TH1D(*aCorrFctn.fRatio);
  66   fkTMonitor = new TH1D(*aCorrFctn.fkTMonitor);
  67 }
  68 //____________________________
  69 AliFemtoQinvCorrFctn::~AliFemtoQinvCorrFctn(){
  70   // destructor
  71   delete fNumerator;
  72   delete fDenominator;
  73   delete fRatio;
  74   delete fkTMonitor;
  75 }
  76 //_________________________
  77 AliFemtoQinvCorrFctn& AliFemtoQinvCorrFctn::operator=(const AliFemtoQinvCorrFctn& aCorrFctn)
  78 {
  79   // assignment operator
  80   if (this == &aCorrFctn)
  81     return *this;
  82
  83   if (fNumerator) delete fNumerator;
  84   fNumerator = new TH1D(*aCorrFctn.fNumerator);
  85   if (fDenominator) delete fDenominator;
  86   fDenominator = new TH1D(*aCorrFctn.fDenominator);
  87   if (fRatio) delete fRatio;
  88   fRatio = new TH1D(*aCorrFctn.fRatio);
  89   if (fkTMonitor) delete fkTMonitor;
  90   fkTMonitor = new TH1D(*aCorrFctn.fkTMonitor);
  91
  92   return *this;
  93 }
  94
  95 //_________________________
  96 void AliFemtoQinvCorrFctn::Finish(){
  97   // here is where we should normalize, fit, etc...
  98   // we should NOT Draw() the histos (as I had done it below),
  99   // since we want to insulate ourselves from root at this level
 100   // of the code.  Do it instead at root command line with browser.
 101   //  fNumerator->Draw();
 102   //fDenominator->Draw();
 103   //fRatio->Draw();
 104   fRatio->Divide(fNumerator,fDenominator,1.0,1.0);
 105
 106 }
 107
 108 //____________________________
 109 AliFemtoString AliFemtoQinvCorrFctn::Report(){
 110   // construct report
 111   string stemp = "Qinv Correlation Function Report:\n";
 112   char ctemp[100];
 113   snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
 114   stemp += ctemp;
 115   snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
 116   stemp += ctemp;
 117   snprintf(ctemp , 100, "Number of entries in ratio:\t%E\n",fRatio->GetEntries());
 118   stemp += ctemp;
 119   //  stemp += mCoulombWeight->Report();
 120   AliFemtoString returnThis = stemp;
 121   return returnThis;
 122 }
 123 //____________________________
 124 void AliFemtoQinvCorrFctn::AddRealPair(AliFemtoPair* pair){
 125   // add true pair
 126   if (fPairCut)
 127     if (!fPairCut->Pass(pair)) return;
 128
 129   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 130   fNumerator->Fill(tQinv);
 131   fkTMonitor->Fill(pair->KT());
 132   //  cout << "AliFemtoQinvCorrFctn::AddRealPair : " << pair->qInv() << " " << tQinv <<
 133   //" " << pair->track1().FourMomentum() << " " << pair->track2().FourMomentum() << endl;
 134 }
 135 //____________________________
 136 void AliFemtoQinvCorrFctn::AddMixedPair(AliFemtoPair* pair){
 137   // add mixed (background) pair
 138   if (fPairCut)
 139     if (!fPairCut->Pass(pair)) return;
 140
 141   double weight = 1.0;
 142   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 143   fDenominator->Fill(tQinv,weight);
 144 }
 145 //____________________________
 146 void AliFemtoQinvCorrFctn::Write(){
 147   // Write out neccessary objects
 148   fNumerator->Write();
 149   fDenominator->Write();
 150   fkTMonitor->Write();
 151 }
 152 //______________________________
 153 TList* AliFemtoQinvCorrFctn::GetOutputList()
 154 {
 155   // Prepare the list of objects to be written to the output
 156   TList *tOutputList = new TList();
 157
 158   tOutputList->Add(fNumerator);
 159   tOutputList->Add(fDenominator);
 160   tOutputList->Add(fkTMonitor);
 161
 162   return tOutputList;
 163 }
 164
 165