PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnNonIdDR.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //                                                                            //
   3 // AliFemtoModelCorrFctnNonIdDR - correlation function for non-identical      //
   4 // particles which uses k* as a function variable. Stores the correlation     //
   5 // function separately for positive and negative signs of k* projections into //
   6 // out, side and long directions, enabling the calculations of double ratios  //
   7 // Uses pair weight to simulate the model correlation function.               //
   8 //                                                                            //
   9 ////////////////////////////////////////////////////////////////////////////////
  10
  11 #include "AliFemtoModelCorrFctnNonIdDR.h"
  12 #include "AliFemtoModelManager.h"
  13 //#include "AliFemtoHisto.h"
  14 #include <cstdio>
  15
  16 #ifdef __ROOT__
  17 ClassImp(AliFemtoModelCorrFctnNonIdDR)
  18 #endif
  19
  20 //____________________________
  21 AliFemtoModelCorrFctnNonIdDR::AliFemtoModelCorrFctnNonIdDR(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
  22   AliFemtoModelCorrFctn(title, nbins, QinvLo, QinvHi),
  23   fNumTOutP(0),
  24   fNumTOutN(0),
  25   fNumTSideP(0),
  26   fNumTSideN(0),
  27   fNumTLongP(0),
  28   fNumTLongN(0),
  29   fNumFOutP(0),
  30   fNumFOutN(0),
  31   fNumFSideP(0),
  32   fNumFSideN(0),
  33   fNumFLongP(0),
  34   fNumFLongN(0),
  35   fDenOutP(0),
  36   fDenOutN(0),
  37   fDenSideP(0),
  38   fDenSideN(0),
  39   fDenLongP(0),
  40   fDenLongN(0)
  41 {
  42   // Default constructor
  43   char bufname[200];
  44
  45   // set up true numerators
  46   snprintf(bufname, 200, "NumTOutP%s", title);
  47   fNumTOutP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  48   snprintf(bufname, 200, "NumTOutN%s", title);
  49   fNumTOutN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  50   snprintf(bufname, 200, "NumTSideP%s", title);
  51   fNumTSideP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  52   snprintf(bufname, 200, "NumTSideN%s", title);
  53   fNumTSideN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  54   snprintf(bufname, 200, "NumTLongP%s", title);
  55   fNumTLongP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  56   snprintf(bufname, 200, "NumTLongN%s", title);
  57   fNumTLongN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  58
  59   // set up fake numerators
  60   snprintf(bufname, 200, "NumFOutP%s", title);
  61   fNumFOutP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  62   snprintf(bufname, 200, "NumFOutN%s", title);
  63   fNumFOutN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  64   snprintf(bufname, 200, "NumFSideP%s", title);
  65   fNumFSideP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  66   snprintf(bufname, 200, "NumFSideN%s", title);
  67   fNumFSideN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  68   snprintf(bufname, 200, "NumFLongP%s", title);
  69   fNumFLongP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  70   snprintf(bufname, 200, "NumFLongN%s", title);
  71   fNumFLongN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  72
  73   // set up denominators
  74   snprintf(bufname, 200, "DenOutP%s", title);
  75   fDenOutP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  76   snprintf(bufname, 200, "DenOutN%s", title);
  77   fDenOutN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  78   snprintf(bufname, 200, "DenSideP%s", title);
  79   fDenSideP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  80   snprintf(bufname, 200, "DenSideN%s", title);
  81   fDenSideN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  82   snprintf(bufname, 200, "DenLongP%s", title);
  83   fDenLongP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  84   snprintf(bufname, 200, "DenLongN%s", title);
  85   fDenLongN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  86
  87   // to enable error bar calculation...
  88   fNumTOutP->Sumw2();
  89   fNumTOutN->Sumw2();
  90   fNumTSideP->Sumw2();
  91   fNumTSideN->Sumw2();
  92   fNumTLongP->Sumw2();
  93   fNumTLongN->Sumw2();
  94   fNumFOutP->Sumw2();
  95   fNumFOutN->Sumw2();
  96   fNumFSideP->Sumw2();
  97   fNumFSideN->Sumw2();
  98   fNumFLongP->Sumw2();
  99   fNumFLongN->Sumw2();
 100   fDenOutP->Sumw2();
 101   fDenOutN->Sumw2();
 102   fDenSideP->Sumw2();
 103   fDenSideN->Sumw2();
 104   fDenLongP->Sumw2();
 105   fDenLongN->Sumw2();
 106 }
 107
 108 //____________________________
 109 AliFemtoModelCorrFctnNonIdDR::AliFemtoModelCorrFctnNonIdDR(const AliFemtoModelCorrFctnNonIdDR& aCorrFctn) :
 110   AliFemtoModelCorrFctn(),
 111   fNumTOutP(0),
 112   fNumTOutN(0),
 113   fNumTSideP(0),
 114   fNumTSideN(0),
 115   fNumTLongP(0),
 116   fNumTLongN(0),
 117   fNumFOutP(0),
 118   fNumFOutN(0),
 119   fNumFSideP(0),
 120   fNumFSideN(0),
 121   fNumFLongP(0),
 122   fNumFLongN(0),
 123   fDenOutP(0),
 124   fDenOutN(0),
 125   fDenSideP(0),
 126   fDenSideN(0),
 127   fDenLongP(0),
 128   fDenLongN(0)
 129 {
 130   // copy constructor
 131   if (aCorrFctn.fNumTOutP)
 132     fNumTOutP = new TH1D(*aCorrFctn.fNumTOutP);
 133   if (aCorrFctn.fNumTOutN)
 134     fNumTOutN = new TH1D(*aCorrFctn.fNumTOutN);
 135   if (aCorrFctn.fNumTSideP)
 136     fNumTSideP = new TH1D(*aCorrFctn.fNumTSideP);
 137   if (aCorrFctn.fNumTSideN)
 138     fNumTSideN = new TH1D(*aCorrFctn.fNumTSideN);
 139   if (aCorrFctn.fNumTLongP)
 140     fNumTLongP = new TH1D(*aCorrFctn.fNumTLongP);
 141   if (aCorrFctn.fNumTLongN)
 142     fNumTLongN = new TH1D(*aCorrFctn.fNumTLongN);
 143
 144   if (aCorrFctn.fNumFOutP)
 145     fNumFOutP = new TH1D(*aCorrFctn.fNumFOutP);
 146   if (aCorrFctn.fNumFOutN)
 147     fNumFOutN = new TH1D(*aCorrFctn.fNumFOutN);
 148   if (aCorrFctn.fNumFSideP)
 149     fNumFSideP = new TH1D(*aCorrFctn.fNumFSideP);
 150   if (aCorrFctn.fNumFSideN)
 151     fNumFSideN = new TH1D(*aCorrFctn.fNumFSideN);
 152   if (aCorrFctn.fNumFLongP)
 153     fNumFLongP = new TH1D(*aCorrFctn.fNumFLongP);
 154   if (aCorrFctn.fNumFLongN)
 155     fNumFLongN = new TH1D(*aCorrFctn.fNumFLongN);
 156
 157   if (aCorrFctn.fDenOutP)
 158     fDenOutP = new TH1D(*aCorrFctn.fDenOutP);
 159   if (aCorrFctn.fDenOutN)
 160     fDenOutN = new TH1D(*aCorrFctn.fDenOutN);
 161   if (aCorrFctn.fDenSideP)
 162     fDenSideP = new TH1D(*aCorrFctn.fDenSideP);
 163   if (aCorrFctn.fDenSideN)
 164     fDenSideN = new TH1D(*aCorrFctn.fDenSideN);
 165   if (aCorrFctn.fDenLongP)
 166     fDenLongP = new TH1D(*aCorrFctn.fDenLongP);
 167   if (aCorrFctn.fDenLongN)
 168     fDenLongN = new TH1D(*aCorrFctn.fDenLongN);
 169 }
 170 //____________________________
 171 AliFemtoModelCorrFctnNonIdDR::~AliFemtoModelCorrFctnNonIdDR(){
 172   // Destructor
 173   delete fNumTOutP;
 174   delete fNumTOutN;
 175   delete fNumTSideP;
 176   delete fNumTSideN;
 177   delete fNumTLongP;
 178   delete fNumTLongN;
 179   delete fNumFOutP;
 180   delete fNumFOutN;
 181   delete fNumFSideP;
 182   delete fNumFSideN;
 183   delete fNumFLongP;
 184   delete fNumFLongN;
 185   delete fDenOutP;
 186   delete fDenOutN;
 187   delete fDenSideP;
 188   delete fDenSideN;
 189   delete fDenLongP;
 190   delete fDenLongN;
 191 }
 192 //_________________________
 193 AliFemtoModelCorrFctnNonIdDR& AliFemtoModelCorrFctnNonIdDR::operator=(const AliFemtoModelCorrFctnNonIdDR& aCorrFctn)
 194 {
 195   // assignment operator
 196   if (this == &aCorrFctn)
 197     return *this;
 198
 199   if (aCorrFctn.fNumTOutP)
 200     fNumTOutP = new TH1D(*aCorrFctn.fNumTOutP);
 201   if (aCorrFctn.fNumTOutN)
 202     fNumTOutN = new TH1D(*aCorrFctn.fNumTOutN);
 203   if (aCorrFctn.fNumTSideP)
 204     fNumTSideP = new TH1D(*aCorrFctn.fNumTSideP);
 205   if (aCorrFctn.fNumTSideN)
 206     fNumTSideN = new TH1D(*aCorrFctn.fNumTSideN);
 207   if (aCorrFctn.fNumTLongP)
 208     fNumTLongP = new TH1D(*aCorrFctn.fNumTLongP);
 209   if (aCorrFctn.fNumTLongN)
 210     fNumTLongN = new TH1D(*aCorrFctn.fNumTLongN);
 211
 212   if (aCorrFctn.fNumFOutP)
 213     fNumFOutP = new TH1D(*aCorrFctn.fNumFOutP);
 214   if (aCorrFctn.fNumFOutN)
 215     fNumFOutN = new TH1D(*aCorrFctn.fNumFOutN);
 216   if (aCorrFctn.fNumFSideP)
 217     fNumFSideP = new TH1D(*aCorrFctn.fNumFSideP);
 218   if (aCorrFctn.fNumFSideN)
 219     fNumFSideN = new TH1D(*aCorrFctn.fNumFSideN);
 220   if (aCorrFctn.fNumFLongP)
 221     fNumFLongP = new TH1D(*aCorrFctn.fNumFLongP);
 222   if (aCorrFctn.fNumFLongN)
 223     fNumFLongN = new TH1D(*aCorrFctn.fNumFLongN);
 224
 225   if (aCorrFctn.fDenOutP)
 226     fDenOutP = new TH1D(*aCorrFctn.fDenOutP);
 227   if (aCorrFctn.fDenOutN)
 228     fDenOutN = new TH1D(*aCorrFctn.fDenOutN);
 229   if (aCorrFctn.fDenSideP)
 230     fDenSideP = new TH1D(*aCorrFctn.fDenSideP);
 231   if (aCorrFctn.fDenSideN)
 232     fDenSideN = new TH1D(*aCorrFctn.fDenSideN);
 233   if (aCorrFctn.fDenLongP)
 234     fDenLongP = new TH1D(*aCorrFctn.fDenLongP);
 235   if (aCorrFctn.fDenLongN)
 236     fDenLongN = new TH1D(*aCorrFctn.fDenLongN);
 237
 238   return *this;
 239 }
 240
 241 //_________________________
 242 void AliFemtoModelCorrFctnNonIdDR::Finish(){
 243   // here is where we should normalize, fit, etc...
 244   // we should NOT Draw() the histos (as I had done it below),
 245   // since we want to insulate ourselves from root at this level
 246   // of the code.  Do it instead at root command line with browser.
 247   //  fNumerator->Draw();
 248   //fDenominator->Draw();
 249   //fRatio->Draw();
 250   //  fRatio->Divide(fNumerator,fDenominator,1.0,1.0);
 251
 252 }
 253
 254 //____________________________
 255 AliFemtoString AliFemtoModelCorrFctnNonIdDR::Report(){
 256   // construct report
 257   string stemp = "Non-identical particles Model Correlation Function Report:\n";
 258   char ctemp[100];
 259   snprintf(ctemp , 100, "Number of entries in numerators:\t%E\n",fNumTOutP->GetEntries()+fNumTOutN->GetEntries());
 260   stemp += ctemp;
 261   snprintf(ctemp , 100, "Number of entries in denominators:\t%E\n",fDenOutP->GetEntries()+fDenOutN->GetEntries());
 262   stemp += ctemp;
 263   //  stemp += mCoulombWeight->Report();
 264   AliFemtoString returnThis = stemp;
 265   return returnThis;
 266 }
 267 //____________________________
 268 void AliFemtoModelCorrFctnNonIdDR::AddRealPair(AliFemtoPair* pair){
 269   // add true pair
 270   double tKStar = pair->KStar();
 271   Double_t weight = fManager->GetWeight(pair);
 272
 273   if (pair->KOut()>0.0)
 274     fNumTOutP->Fill(tKStar, weight);
 275   else
 276     fNumTOutN->Fill(tKStar, weight);
 277
 278   if (pair->KSide()>0.0)
 279     fNumTSideP->Fill(tKStar, weight);
 280   else
 281     fNumTSideN->Fill(tKStar, weight);
 282
 283   if (pair->KLong()>0.0)
 284     fNumTLongP->Fill(tKStar, weight);
 285   else
 286     fNumTLongN->Fill(tKStar, weight);
 287
 288 }
 289 //____________________________
 290 void AliFemtoModelCorrFctnNonIdDR::AddMixedPair(AliFemtoPair* pair){
 291   // add mixed (background) pair
 292   double tKStar = pair->KStar();
 293   Double_t weight = fManager->GetWeight(pair);
 294
 295   if (pair->KOut()>0.0)
 296     fNumFOutP->Fill(tKStar, weight);
 297   else
 298     fNumFOutN->Fill(tKStar, weight);
 299
 300   if (pair->KSide()>0.0)
 301     fNumFSideP->Fill(tKStar, weight);
 302   else
 303     fNumFSideN->Fill(tKStar, weight);
 304
 305   if (pair->KLong()>0.0)
 306     fNumFLongP->Fill(tKStar, weight);
 307   else
 308     fNumFLongN->Fill(tKStar, weight);
 309
 310   if (pair->KOut()>0.0)
 311     fDenOutP->Fill(tKStar);
 312   else
 313     fDenOutN->Fill(tKStar);
 314
 315   if (pair->KSide()>0.0)
 316     fDenSideP->Fill(tKStar);
 317   else
 318     fDenSideN->Fill(tKStar);
 319
 320   if (pair->KLong()>0.0)
 321     fDenLongP->Fill(tKStar);
 322   else
 323     fDenLongN->Fill(tKStar);
 324 }
 325 //____________________________
 326 void AliFemtoModelCorrFctnNonIdDR::Write(){
 327   // Write out histos
 328   fNumTOutP->Write();
 329   fNumTOutN->Write();
 330   fNumTSideP->Write();
 331   fNumTSideN->Write();
 332   fNumTLongP->Write();
 333   fNumTLongN->Write();
 334   fNumFOutP->Write();
 335   fNumFOutN->Write();
 336   fNumFSideP->Write();
 337   fNumFSideN->Write();
 338   fNumFLongP->Write();
 339   fNumFLongN->Write();
 340   fDenOutP->Write();
 341   fDenOutN->Write();
 342   fDenSideP->Write();
 343   fDenSideN->Write();
 344   fDenLongP->Write();
 345   fDenLongN->Write();
 346 }
 347
 348 TList* AliFemtoModelCorrFctnNonIdDR::GetOutputList()
 349 {
 350   // Prepare the list of objects to be written to the output
 351   TList *tOutputList = new TList();
 352
 353   tOutputList->Add(fNumTOutP);
 354   tOutputList->Add(fNumTOutN);
 355   tOutputList->Add(fNumTSideP);
 356   tOutputList->Add(fNumTSideN);
 357   tOutputList->Add(fNumTLongP);
 358   tOutputList->Add(fNumTLongN);
 359   tOutputList->Add(fNumFOutP);
 360   tOutputList->Add(fNumFOutN);
 361   tOutputList->Add(fNumFSideP);
 362   tOutputList->Add(fNumFSideN);
 363   tOutputList->Add(fNumFLongP);
 364   tOutputList->Add(fNumFLongN);
 365   tOutputList->Add(fDenOutP);
 366   tOutputList->Add(fDenOutN);
 367   tOutputList->Add(fDenSideP);
 368   tOutputList->Add(fDenSideN);
 369   tOutputList->Add(fDenLongP);
 370   tOutputList->Add(fDenLongN);
 371
 372   return tOutputList;
 373 }
 374
 375 //_______________________
 376 AliFemtoModelCorrFctn* AliFemtoModelCorrFctnNonIdDR::Clone()
 377 {
 378   // Create clone
 379   AliFemtoModelCorrFctnNonIdDR *tCopy = new AliFemtoModelCorrFctnNonIdDR(*this);
 380
 381   return tCopy;
 382 }