PWGCF/FEMTOSCOPY/AliFemto/AliFemtoCorrFctnNonIdDR.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //                                                                            //
   3 // AliFemtoCorrFctnNonIdDR - correlation function for non-identical particles //
   4 // uses k* as a function variable. Stores the correlation function separately //
   5 // for positive and negative signs of k* projections into out, side and long  //
   6 // directions, enabling the calculations of double ratios                     //
   7 //                                                                            //
   8 ////////////////////////////////////////////////////////////////////////////////
   9
  10 #include "AliFemtoCorrFctnNonIdDR.h"
  11 //#include "AliFemtoHisto.h"
  12 #include <cstdio>
  13
  14 #ifdef __ROOT__
  15 ClassImp(AliFemtoCorrFctnNonIdDR)
  16 #endif
  17
  18 //____________________________
  19 AliFemtoCorrFctnNonIdDR::AliFemtoCorrFctnNonIdDR(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
  20   fNumOutP(0),
  21   fNumOutN(0),
  22   fNumSideP(0),
  23   fNumSideN(0),
  24   fNumLongP(0),
  25   fNumLongN(0),
  26   fDenOutP(0),
  27   fDenOutN(0),
  28   fDenSideP(0),
  29   fDenSideN(0),
  30   fDenLongP(0),
  31     fDenLongN(0),
  32     fkTMonitor(0)
  33
  34 {
  35   // Default constructor
  36   // set up numerators
  37   char bufname[200];
  38   snprintf(bufname, 200, "NumOutP%s", title);
  39   fNumOutP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  40   snprintf(bufname, 200, "NumOutN%s", title);
  41   fNumOutN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  42   snprintf(bufname, 200, "NumSideP%s", title);
  43   fNumSideP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  44   snprintf(bufname, 200, "NumSideN%s", title);
  45   fNumSideN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  46   snprintf(bufname, 200, "NumLongP%s", title);
  47   fNumLongP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  48   snprintf(bufname, 200, "NumLongN%s", title);
  49   fNumLongN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  50
  51   // set up denominators
  52   snprintf(bufname, 200, "DenOutP%s", title);
  53   fDenOutP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  54   snprintf(bufname, 200, "DenOutN%s", title);
  55   fDenOutN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  56   snprintf(bufname, 200, "DenSideP%s", title);
  57   fDenSideP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  58   snprintf(bufname, 200, "DenSideN%s", title);
  59   fDenSideN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  60   snprintf(bufname, 200, "DenLongP%s", title);
  61   fDenLongP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  62   snprintf(bufname, 200, "DenLongN%s", title);
  63   fDenLongN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  64
  65     char tTitkT[101] = "kTDep";
  66     strncat(tTitkT,title, 100);
  67     fkTMonitor = new TH1D(tTitkT,title,250,0.0,5.0);
  68
  69   // to enable error bar calculation...
  70   fNumOutP->Sumw2();
  71   fNumOutN->Sumw2();
  72   fNumSideP->Sumw2();
  73   fNumSideN->Sumw2();
  74   fNumLongP->Sumw2();
  75   fNumLongN->Sumw2();
  76   fDenOutP->Sumw2();
  77   fDenOutN->Sumw2();
  78   fDenSideP->Sumw2();
  79   fDenSideN->Sumw2();
  80   fDenLongP->Sumw2();
  81   fDenLongN->Sumw2();
  82
  83     fkTMonitor->Sumw2();
  84 }
  85
  86 //____________________________
  87 AliFemtoCorrFctnNonIdDR::AliFemtoCorrFctnNonIdDR(const AliFemtoCorrFctnNonIdDR& aCorrFctn) :
  88   AliFemtoCorrFctn(),
  89   fNumOutP(0),
  90   fNumOutN(0),
  91   fNumSideP(0),
  92   fNumSideN(0),
  93   fNumLongP(0),
  94   fNumLongN(0),
  95   fDenOutP(0),
  96   fDenOutN(0),
  97   fDenSideP(0),
  98   fDenSideN(0),
  99   fDenLongP(0),
 100     fDenLongN(0),
 101     fkTMonitor(0)
 102 {
 103   // copy constructor
 104   if (aCorrFctn.fNumOutP)
 105     fNumOutP = new TH1D(*aCorrFctn.fNumOutP);
 106   if (aCorrFctn.fNumOutN)
 107     fNumOutN = new TH1D(*aCorrFctn.fNumOutN);
 108   if (aCorrFctn.fNumSideP)
 109     fNumSideP = new TH1D(*aCorrFctn.fNumSideP);
 110   if (aCorrFctn.fNumSideN)
 111     fNumSideN = new TH1D(*aCorrFctn.fNumSideN);
 112   if (aCorrFctn.fNumLongP)
 113     fNumLongP = new TH1D(*aCorrFctn.fNumLongP);
 114   if (aCorrFctn.fNumLongN)
 115     fNumLongN = new TH1D(*aCorrFctn.fNumLongN);
 116
 117   if (aCorrFctn.fDenOutP)
 118     fDenOutP = new TH1D(*aCorrFctn.fDenOutP);
 119   if (aCorrFctn.fDenOutN)
 120     fDenOutN = new TH1D(*aCorrFctn.fDenOutN);
 121   if (aCorrFctn.fDenSideP)
 122     fDenSideP = new TH1D(*aCorrFctn.fDenSideP);
 123   if (aCorrFctn.fDenSideN)
 124     fDenSideN = new TH1D(*aCorrFctn.fDenSideN);
 125   if (aCorrFctn.fDenLongP)
 126     fDenLongP = new TH1D(*aCorrFctn.fDenLongP);
 127   if (aCorrFctn.fDenLongN)
 128     fDenLongN = new TH1D(*aCorrFctn.fDenLongN);
 129
 130     if (aCorrFctn.fkTMonitor)
 131         fkTMonitor = new TH1D(*aCorrFctn.fkTMonitor);
 132 }
 133
 134 //____________________________
 135 AliFemtoCorrFctnNonIdDR::~AliFemtoCorrFctnNonIdDR(){
 136   delete fNumOutP;
 137   delete fNumOutN;
 138   delete fNumSideP;
 139   delete fNumSideN;
 140   delete fNumLongP;
 141   delete fNumLongN;
 142   delete fDenOutP;
 143   delete fDenOutN;
 144   delete fDenSideP;
 145   delete fDenSideN;
 146   delete fDenLongP;
 147   delete fDenLongN;
 148     delete fkTMonitor;
 149 }
 150
 151 //_________________________
 152 AliFemtoCorrFctnNonIdDR& AliFemtoCorrFctnNonIdDR::operator=(const AliFemtoCorrFctnNonIdDR& aCorrFctn)
 153 {
 154   // assignment operator
 155   if (this == &aCorrFctn)
 156     return *this;
 157
 158   if (aCorrFctn.fNumOutP)
 159     fNumOutP = new TH1D(*aCorrFctn.fNumOutP);
 160   if (aCorrFctn.fNumOutN)
 161     fNumOutN = new TH1D(*aCorrFctn.fNumOutN);
 162   if (aCorrFctn.fNumSideP)
 163     fNumSideP = new TH1D(*aCorrFctn.fNumSideP);
 164   if (aCorrFctn.fNumSideN)
 165     fNumSideN = new TH1D(*aCorrFctn.fNumSideN);
 166   if (aCorrFctn.fNumLongP)
 167     fNumLongP = new TH1D(*aCorrFctn.fNumLongP);
 168   if (aCorrFctn.fNumLongN)
 169     fNumLongN = new TH1D(*aCorrFctn.fNumLongN);
 170
 171   if (aCorrFctn.fDenOutP)
 172     fDenOutP = new TH1D(*aCorrFctn.fDenOutP);
 173   if (aCorrFctn.fDenOutN)
 174     fDenOutN = new TH1D(*aCorrFctn.fDenOutN);
 175   if (aCorrFctn.fDenSideP)
 176     fDenSideP = new TH1D(*aCorrFctn.fDenSideP);
 177   if (aCorrFctn.fDenSideN)
 178     fDenSideN = new TH1D(*aCorrFctn.fDenSideN);
 179   if (aCorrFctn.fDenLongP)
 180     fDenLongP = new TH1D(*aCorrFctn.fDenLongP);
 181   if (aCorrFctn.fDenLongN)
 182     fDenLongN = new TH1D(*aCorrFctn.fDenLongN);
 183
 184     if (aCorrFctn.fkTMonitor)
 185         fkTMonitor = new TH1D(*aCorrFctn.fkTMonitor);
 186
 187   return *this;
 188 }
 189
 190 //_________________________
 191 void AliFemtoCorrFctnNonIdDR::Finish(){
 192   // here is where we should normalize, fit, etc...
 193   // we should NOT Draw() the histos (as I had done it below),
 194   // since we want to insulate ourselves from root at this level
 195   // of the code.  Do it instead at root command line with browser.
 196   //  fNumerator->Draw();
 197   //fDenominator->Draw();
 198   //fRatio->Draw();
 199   //  fRatio->Divide(fNumerator,fDenominator,1.0,1.0);
 200
 201 }
 202
 203 //____________________________
 204 AliFemtoString AliFemtoCorrFctnNonIdDR::Report(){
 205   // construct report
 206   string stemp = "Non-identical particles Correlation Function Report:\n";
 207   char ctemp[100];
 208   snprintf(ctemp , 100, "Number of entries in numerators:\t%E\n",fNumOutP->GetEntries()+fNumOutN->GetEntries());
 209   stemp += ctemp;
 210   snprintf(ctemp , 100, "Number of entries in denominators:\t%E\n",fDenOutP->GetEntries()+fDenOutN->GetEntries());
 211   stemp += ctemp;
 212   //  stemp += mCoulombWeight->Report();
 213   AliFemtoString returnThis = stemp;
 214   return returnThis;
 215 }
 216 //____________________________
 217 void AliFemtoCorrFctnNonIdDR::AddRealPair(AliFemtoPair* pair){
 218   // add true pair
 219
 220     if (fPairCut)
 221         if (!fPairCut->Pass(pair)) return;
 222
 223   double tKStar = pair->KStar();
 224   if (pair->KOut()>0.0)
 225     fNumOutP->Fill(tKStar);
 226   else
 227     fNumOutN->Fill(tKStar);
 228
 229   if (pair->KSide()>0.0)
 230     fNumSideP->Fill(tKStar);
 231   else
 232     fNumSideN->Fill(tKStar);
 233
 234   if (pair->KLong()>0.0)
 235     fNumLongP->Fill(tKStar);
 236   else
 237     fNumLongN->Fill(tKStar);
 238
 239     fkTMonitor->Fill(pair->KT());
 240
 241 }
 242 //____________________________
 243 void AliFemtoCorrFctnNonIdDR::AddMixedPair(AliFemtoPair* pair){
 244   // add mixed (background) pair
 245
 246     if (fPairCut)
 247         if (!fPairCut->Pass(pair)) return;
 248
 249   double tKStar = pair->KStar();
 250   if (pair->KOut()>0.0)
 251     fDenOutP->Fill(tKStar);
 252   else
 253     fDenOutN->Fill(tKStar);
 254
 255   if (pair->KSide()>0.0)
 256     fDenSideP->Fill(tKStar);
 257   else
 258     fDenSideN->Fill(tKStar);
 259
 260   if (pair->KLong()>0.0)
 261     fDenLongP->Fill(tKStar);
 262   else
 263     fDenLongN->Fill(tKStar);
 264 }
 265 //____________________________
 266 void AliFemtoCorrFctnNonIdDR::Write(){
 267   fNumOutP->Write();
 268   fNumOutN->Write();
 269   fNumSideP->Write();
 270   fNumSideN->Write();
 271   fNumLongP->Write();
 272   fNumLongN->Write();
 273   fDenOutP->Write();
 274   fDenOutN->Write();
 275   fDenSideP->Write();
 276   fDenSideN->Write();
 277   fDenLongP->Write();
 278   fDenLongN->Write();
 279     fkTMonitor->Write();
 280 }
 281
 282 TList* AliFemtoCorrFctnNonIdDR::GetOutputList()
 283 {
 284   // Prepare the list of objects to be written to the output
 285   TList *tOutputList = new TList();
 286
 287   tOutputList->Add(fNumOutP);
 288   tOutputList->Add(fNumOutN);
 289   tOutputList->Add(fNumSideP);
 290   tOutputList->Add(fNumSideN);
 291   tOutputList->Add(fNumLongP);
 292   tOutputList->Add(fNumLongN);
 293   tOutputList->Add(fDenOutP);
 294   tOutputList->Add(fDenOutN);
 295   tOutputList->Add(fDenSideP);
 296   tOutputList->Add(fDenSideN);
 297   tOutputList->Add(fDenLongP);
 298   tOutputList->Add(fDenLongN);
 299     tOutputList->Add(fkTMonitor);
 300
 301   return tOutputList;
 302 }