PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnDEtaDPhi.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //                                                                            //
   3 // AliFemtoCorrFctnDEtaDPhi - A correlation function that analyzes            //
   4 // two particle correlations with respect to the azimuthal angle (phi)        //
   5 // and pseudorapidity (eta) difference                                        //
   6 //                                                                            //
   7 // Authors: Adam Kisiel Adam.Kisiel@cern.ch                                   //
   8 //                                                                            //
   9 ////////////////////////////////////////////////////////////////////////////////
  10
  11 #include "AliFemtoCorrFctnDEtaDPhi.h"
  12 //#include "AliFemtoHisto.hh"
  13 #include <cstdio>
  14 #include <TMath.h>
  15
  16 #ifdef __ROOT__
  17 ClassImp(AliFemtoCorrFctnDEtaDPhi)
  18 #endif
  19
  20 //____________________________
  21 AliFemtoCorrFctnDEtaDPhi::AliFemtoCorrFctnDEtaDPhi(char* title, const int& aPhiBins=20, const int& aEtaBins=20):
  22   AliFemtoCorrFctn(),
  23   fDPhiDEtaNumerator(0),
  24   fDPhiDEtaDenominator(0),
  25   fDPhiDEtaColNumerator(0),
  26   fDPhiDEtaColDenominator(0),
  27   fDPhiNumerator(0),
  28   fDPhiDenominator(0),
  29   fDCosNumerator(0),
  30   fDCosDenominator(0),
  31   fDPhiPtNumerator(0),
  32   fDPhiPtDenominator(0),
  33   fDCosPtNumerator(0),
  34   fDCosPtDenominator(0)
  35 {
  36   // set up numerator
  37   char tTitNumD[100] = "NumDPhiDEta";
  38   strcat(tTitNumD,title);
  39   fDPhiDEtaNumerator = new TH2D(tTitNumD,title,aPhiBins,-0.5*TMath::Pi(),1.5*TMath::Pi(),aEtaBins,-2.0,2.0);
  40   // set up denominator
  41   char tTitDenD[100] = "DenDPhiDEta";
  42   strcat(tTitDenD,title);
  43   fDPhiDEtaDenominator = new TH2D(tTitDenD,title,aPhiBins,-0.5*TMath::Pi(),1.5*TMath::Pi(),aEtaBins,-2.0,2.0);
  44
  45   // set up numerator
  46   char tTitNumR[100] = "NumDPhiDEtaCol";
  47   strcat(tTitNumR,title);
  48   fDPhiDEtaColNumerator = new TH2D(tTitNumR,title,aPhiBins,-0.5*TMath::Pi(),1.5*TMath::Pi(),aEtaBins,-2.0,2.0);
  49   // set up denominator
  50   char tTitDenR[100] = "DenDPhiDEtaCol";
  51   strcat(tTitDenR,title);
  52   fDPhiDEtaColDenominator = new TH2D(tTitDenR,title,aPhiBins,-0.5*TMath::Pi(),1.5*TMath::Pi(),aEtaBins,-2.0,2.0);
  53
  54   // set up numerator
  55   char tTitNumDPhi[100] = "NumDPhi";
  56   strcat(tTitNumDPhi,title);
  57   fDPhiNumerator = new TH1D(tTitNumDPhi,title,aPhiBins*2,-0.5*TMath::Pi(),1.5*TMath::Pi());
  58   // set up denominator
  59   char tTitDenDPhi[100] = "DenDPhi";
  60   strcat(tTitDenDPhi,title);
  61   fDPhiDenominator = new TH1D(tTitDenDPhi,title,aPhiBins*2,-0.5*TMath::Pi(),1.5*TMath::Pi());
  62
  63   // set up numerator
  64   char tTitNumDCos[100] = "NumDCos";
  65   strcat(tTitNumDCos,title);
  66   fDCosNumerator = new TH1D(tTitNumDCos,title,aPhiBins*2,-1.0,1.0);
  67   // set up denominator
  68   char tTitDenDCos[100] = "DenDCos";
  69   strcat(tTitDenDCos,title);
  70   fDCosDenominator = new TH1D(tTitDenDCos,title,aPhiBins*2,-1.0,1.0);
  71
  72   // set up numerator
  73   char tTitNumDPhiPt[100] = "NumDPhiPt";
  74   strcat(tTitNumDPhiPt,title);
  75   fDPhiPtNumerator = new TH2D(tTitNumDPhiPt,title,aPhiBins*2,-0.5*TMath::Pi(),1.5*TMath::Pi(), 30, 0.0, 3.0);
  76   // set up denominator
  77   char tTitDenDPhiPt[100] = "DenDPhiPt";
  78   strcat(tTitDenDPhiPt,title);
  79   fDPhiPtDenominator = new TH2D(tTitDenDPhiPt,title,aPhiBins*2,-0.5*TMath::Pi(),1.5*TMath::Pi(), 30, 0.0, 3.0);
  80
  81   // set up numerator
  82   char tTitNumDCosPt[100] = "NumDCosPt";
  83   strcat(tTitNumDCosPt,title);
  84   fDCosPtNumerator = new TH2D(tTitNumDCosPt,title,aPhiBins*2,-1.0,1.0, 30, 0.0, 3.0);
  85   // set up denominator
  86   char tTitDenDCosPt[100] = "DenDCosPt";
  87   strcat(tTitDenDCosPt,title);
  88   fDCosPtDenominator = new TH2D(tTitDenDCosPt,title,aPhiBins*2,-1.0,1.0, 30, 0.0, 3.0);
  89
  90   // to enable error bar calculation...
  91   fDPhiDEtaNumerator->Sumw2();
  92   fDPhiDEtaDenominator->Sumw2();
  93   fDPhiDEtaColNumerator->Sumw2();
  94   fDPhiDEtaColDenominator->Sumw2();
  95   fDPhiNumerator->Sumw2();
  96   fDPhiDenominator->Sumw2();
  97   fDCosNumerator->Sumw2();
  98   fDCosDenominator->Sumw2();
  99   fDPhiPtNumerator->Sumw2();
 100   fDPhiPtDenominator->Sumw2();
 101   fDCosPtNumerator->Sumw2();
 102   fDCosPtDenominator->Sumw2();
 103
 104 }
 105
 106 //____________________________
 107 AliFemtoCorrFctnDEtaDPhi::AliFemtoCorrFctnDEtaDPhi(const AliFemtoCorrFctnDEtaDPhi& aCorrFctn) :
 108   AliFemtoCorrFctn(),
 109   fDPhiDEtaNumerator(0),
 110   fDPhiDEtaDenominator(0),
 111   fDPhiDEtaColNumerator(0),
 112   fDPhiDEtaColDenominator(0),
 113   fDPhiNumerator(0),
 114   fDPhiDenominator(0),
 115   fDCosNumerator(0),
 116   fDCosDenominator(0),
 117   fDPhiPtNumerator(0),
 118   fDPhiPtDenominator(0),
 119   fDCosPtNumerator(0),
 120   fDCosPtDenominator(0)
 121 {
 122   // copy constructor
 123   if (aCorrFctn.fDPhiDEtaNumerator)
 124     fDPhiDEtaNumerator = new TH2D(*aCorrFctn.fDPhiDEtaNumerator);
 125   else
 126     fDPhiDEtaNumerator = 0;
 127   if (aCorrFctn.fDPhiDEtaDenominator)
 128     fDPhiDEtaDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 129   else
 130     fDPhiDEtaDenominator = 0;
 131
 132   if (aCorrFctn.fDPhiDEtaColNumerator)
 133     fDPhiDEtaColNumerator = new TH2D(*aCorrFctn.fDPhiDEtaColNumerator);
 134   else
 135     fDPhiDEtaColNumerator = 0;
 136   if (aCorrFctn.fDPhiDEtaColDenominator)
 137     fDPhiDEtaColDenominator = new TH2D(*aCorrFctn.fDPhiDEtaColDenominator);
 138   else
 139     fDPhiDEtaColDenominator = 0;
 140
 141   if (aCorrFctn.fDPhiNumerator)
 142     fDPhiNumerator = new TH1D(*aCorrFctn.fDPhiNumerator);
 143   else
 144     fDPhiNumerator = 0;
 145   if (aCorrFctn.fDPhiDenominator)
 146     fDPhiDenominator = new TH1D(*aCorrFctn.fDPhiDenominator);
 147   else
 148     fDPhiDenominator = 0;
 149
 150   if (aCorrFctn.fDCosNumerator)
 151     fDCosNumerator = new TH1D(*aCorrFctn.fDCosNumerator);
 152   else
 153     fDCosNumerator = 0;
 154   if (aCorrFctn.fDCosDenominator)
 155     fDCosDenominator = new TH1D(*aCorrFctn.fDCosDenominator);
 156   else
 157     fDCosDenominator = 0;
 158
 159   if (aCorrFctn.fDPhiPtNumerator)
 160     fDPhiPtNumerator = new TH2D(*aCorrFctn.fDPhiPtNumerator);
 161   else
 162     fDPhiPtNumerator = 0;
 163   if (aCorrFctn.fDPhiPtDenominator)
 164     fDPhiPtDenominator = new TH2D(*aCorrFctn.fDPhiPtDenominator);
 165   else
 166     fDPhiPtDenominator = 0;
 167
 168   if (aCorrFctn.fDCosPtNumerator)
 169     fDCosPtNumerator = new TH2D(*aCorrFctn.fDCosPtNumerator);
 170   else
 171     fDCosPtNumerator = 0;
 172   if (aCorrFctn.fDCosPtDenominator)
 173     fDCosPtDenominator = new TH2D(*aCorrFctn.fDCosPtDenominator);
 174   else
 175     fDCosPtDenominator = 0;
 176 }
 177 //____________________________
 178 AliFemtoCorrFctnDEtaDPhi::~AliFemtoCorrFctnDEtaDPhi(){
 179   // destructor
 180   delete fDPhiDEtaNumerator;
 181   delete fDPhiDEtaDenominator;
 182   delete fDPhiDEtaColNumerator;
 183   delete fDPhiDEtaColDenominator;
 184   delete fDPhiNumerator;
 185   delete fDPhiDenominator;
 186   delete fDCosNumerator;
 187   delete fDCosDenominator;
 188   delete fDPhiPtNumerator;
 189   delete fDPhiPtDenominator;
 190   delete fDCosPtNumerator;
 191   delete fDCosPtDenominator;
 192 }
 193 //_________________________
 194 AliFemtoCorrFctnDEtaDPhi& AliFemtoCorrFctnDEtaDPhi::operator=(const AliFemtoCorrFctnDEtaDPhi& aCorrFctn)
 195 {
 196   // assignment operator
 197   if (this == &aCorrFctn)
 198     return *this;
 199
 200   if (aCorrFctn.fDPhiDEtaNumerator)
 201     fDPhiDEtaNumerator = new TH2D(*aCorrFctn.fDPhiDEtaNumerator);
 202   else
 203     fDPhiDEtaNumerator = 0;
 204   if (aCorrFctn.fDPhiDEtaDenominator)
 205     fDPhiDEtaDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 206   else
 207     fDPhiDEtaDenominator = 0;
 208
 209   if (aCorrFctn.fDPhiDEtaColNumerator)
 210     fDPhiDEtaColNumerator = new TH2D(*aCorrFctn.fDPhiDEtaColNumerator);
 211   else
 212     fDPhiDEtaColNumerator = 0;
 213   if (aCorrFctn.fDPhiDEtaColDenominator)
 214     fDPhiDEtaColDenominator = new TH2D(*aCorrFctn.fDPhiDEtaColDenominator);
 215   else
 216     fDPhiDEtaColDenominator = 0;
 217
 218   if (aCorrFctn.fDPhiNumerator)
 219     fDPhiNumerator = new TH1D(*aCorrFctn.fDPhiNumerator);
 220   else
 221     fDPhiNumerator = 0;
 222   if (aCorrFctn.fDPhiDenominator)
 223     fDPhiDenominator = new TH1D(*aCorrFctn.fDPhiDenominator);
 224   else
 225     fDPhiDenominator = 0;
 226
 227   if (aCorrFctn.fDCosNumerator)
 228     fDCosNumerator = new TH1D(*aCorrFctn.fDCosNumerator);
 229   else
 230     fDCosNumerator = 0;
 231   if (aCorrFctn.fDCosDenominator)
 232     fDCosDenominator = new TH1D(*aCorrFctn.fDCosDenominator);
 233   else
 234     fDCosDenominator = 0;
 235
 236   if (aCorrFctn.fDPhiPtNumerator)
 237     fDPhiPtNumerator = new TH2D(*aCorrFctn.fDPhiPtNumerator);
 238   else
 239     fDPhiPtNumerator = 0;
 240   if (aCorrFctn.fDPhiPtDenominator)
 241     fDPhiPtDenominator = new TH2D(*aCorrFctn.fDPhiPtDenominator);
 242   else
 243     fDPhiPtDenominator = 0;
 244
 245   if (aCorrFctn.fDCosPtNumerator)
 246     fDCosPtNumerator = new TH2D(*aCorrFctn.fDCosPtNumerator);
 247   else
 248     fDCosPtNumerator = 0;
 249   if (aCorrFctn.fDCosPtDenominator)
 250     fDCosPtDenominator = new TH2D(*aCorrFctn.fDCosPtDenominator);
 251   else
 252     fDCosPtDenominator = 0;
 253
 254   return *this;
 255 }
 256 //_________________________
 257 void AliFemtoCorrFctnDEtaDPhi::Finish(){
 258   // here is where we should normalize, fit, etc...
 259   // we should NOT Draw() the histos (as I had done it below),
 260   // since we want to insulate ourselves from root at this level
 261   // of the code.  Do it instead at root command line with browser.
 262   //  mShareNumerator->Draw();
 263   //mShareDenominator->Draw();
 264   //mRatio->Draw();
 265
 266 }
 267
 268 //____________________________
 269 AliFemtoString AliFemtoCorrFctnDEtaDPhi::Report(){
 270   // create report
 271   string stemp = "TPC Ncls Correlation Function Report:\n";
 272   char ctemp[100];
 273   sprintf(ctemp,"Number of entries in numerator:\t%E\n",fDPhiDEtaNumerator->GetEntries());
 274   stemp += ctemp;
 275   sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDPhiDEtaDenominator->GetEntries());
 276   stemp += ctemp;
 277   //  stemp += mCoulombWeight->Report();
 278   AliFemtoString returnThis = stemp;
 279   return returnThis;
 280 }
 281 //____________________________
 282 void AliFemtoCorrFctnDEtaDPhi::AddRealPair( AliFemtoPair* pair){
 283   // add real (effect) pair
 284   double phi1 = pair->Track1()->Track()->P().phi();
 285   double phi2 = pair->Track2()->Track()->P().phi();
 286   double eta1 = pair->Track1()->Track()->P().pseudoRapidity();
 287   double eta2 = pair->Track2()->Track()->P().pseudoRapidity();
 288
 289   double dphi = phi1 - phi2;
 290   while (dphi<-TMath::Pi()/2) dphi+=TMath::Pi()*2;
 291   while (dphi>3*TMath::Pi()/2) dphi-=TMath::Pi()*2;
 292
 293   double deta = eta1 - eta2;
 294
 295   double px1 = pair->Track1()->Track()->P().x();
 296   double py1 = pair->Track1()->Track()->P().y();
 297   double pz1 = pair->Track1()->Track()->P().z();
 298
 299   double px2 = pair->Track2()->Track()->P().x();
 300   double py2 = pair->Track2()->Track()->P().y();
 301   double pz2 = pair->Track2()->Track()->P().z();
 302
 303   double pt1 = TMath::Hypot(px1, py1);
 304   double pt2 = TMath::Hypot(px2, py2);
 305   double ptmin = pt1>pt2 ? pt2 : pt1;
 306
 307   double cosphi = (px1*px2 + py1*py2 + pz1*pz2)/
 308     sqrt((px1*px1 + py1*py1 + pz1*pz1)*(px2*px2 + py2*py2 + pz2*pz2));
 309
 310   fDPhiDEtaNumerator->Fill(dphi, deta);
 311
 312   if (cosphi > 0) {
 313     fDPhiDEtaColNumerator->Fill(dphi, deta);
 314   }
 315   else {
 316     fDPhiDEtaColNumerator->Fill(dphi, -eta1-eta2);
 317   }
 318
 319   fDPhiNumerator->Fill(dphi);
 320   fDCosNumerator->Fill(cosphi);
 321
 322   fDPhiPtNumerator->Fill(dphi, ptmin);
 323   fDCosPtNumerator->Fill(cosphi, ptmin);
 324
 325 }
 326 //____________________________
 327 void AliFemtoCorrFctnDEtaDPhi::AddMixedPair( AliFemtoPair* pair){
 328   // add mixed (background) pair
 329   double phi1 = pair->Track1()->Track()->P().phi();
 330   double phi2 = pair->Track2()->Track()->P().phi();
 331   double eta1 = pair->Track1()->Track()->P().pseudoRapidity();
 332   double eta2 = pair->Track2()->Track()->P().pseudoRapidity();
 333
 334   double dphi = phi1 - phi2;
 335   while (dphi<-TMath::Pi()/2) dphi+=TMath::Pi()*2;
 336   while (dphi>3*TMath::Pi()/2) dphi-=TMath::Pi()*2;
 337
 338   double deta = eta1 - eta2;
 339
 340   double px1 = pair->Track1()->Track()->P().x();
 341   double py1 = pair->Track1()->Track()->P().y();
 342   double pz1 = pair->Track1()->Track()->P().z();
 343
 344   double px2 = pair->Track2()->Track()->P().x();
 345   double py2 = pair->Track2()->Track()->P().y();
 346   double pz2 = pair->Track2()->Track()->P().z();
 347
 348   double pt1 = TMath::Hypot(px1, py1);
 349   double pt2 = TMath::Hypot(px2, py2);
 350   double ptmin = pt1>pt2 ? pt2 : pt1;
 351
 352   double cosphi = (px1*px2 + py1*py2 + pz1*pz2)/
 353     sqrt((px1*px1 + py1*py1 + pz1*pz1)*(px2*px2 + py2*py2 + pz2*pz2));
 354
 355   fDPhiDEtaDenominator->Fill(dphi, deta);
 356
 357   if (cosphi > 0) {
 358     fDPhiDEtaColDenominator->Fill(dphi, deta);
 359   }
 360   else {
 361     fDPhiDEtaColDenominator->Fill(dphi, -eta1-eta2);
 362   }
 363
 364   fDPhiDenominator->Fill(dphi);
 365   fDCosDenominator->Fill(cosphi);
 366
 367   fDPhiPtDenominator->Fill(dphi, ptmin);
 368   fDCosPtDenominator->Fill(cosphi, ptmin);
 369 }
 370
 371
 372 void AliFemtoCorrFctnDEtaDPhi::WriteHistos()
 373 {
 374   // Write out result histograms
 375   fDPhiDEtaNumerator->Write();
 376   fDPhiDEtaDenominator->Write();
 377   fDPhiDEtaColNumerator->Write();
 378   fDPhiDEtaColDenominator->Write();
 379   fDPhiNumerator->Write();
 380   fDPhiDenominator->Write();
 381   fDCosNumerator->Write();
 382   fDCosDenominator->Write();
 383   fDPhiPtNumerator->Write();
 384   fDPhiPtDenominator->Write();
 385   fDCosPtNumerator->Write();
 386   fDCosPtDenominator->Write();
 387 }
 388
 389 TList* AliFemtoCorrFctnDEtaDPhi::GetOutputList()
 390 {
 391   // Prepare the list of objects to be written to the output
 392   TList *tOutputList = new TList();
 393
 394   tOutputList->Add(fDPhiDEtaNumerator);
 395   tOutputList->Add(fDPhiDEtaDenominator);
 396   tOutputList->Add(fDPhiDEtaColNumerator);
 397   tOutputList->Add(fDPhiDEtaColDenominator);
 398   tOutputList->Add(fDPhiNumerator);
 399   tOutputList->Add(fDPhiDenominator);
 400   tOutputList->Add(fDCosNumerator);
 401   tOutputList->Add(fDCosDenominator);
 402   tOutputList->Add(fDPhiPtNumerator);
 403   tOutputList->Add(fDPhiPtDenominator);
 404   tOutputList->Add(fDCosPtNumerator);
 405   tOutputList->Add(fDCosPtDenominator);
 406
 407   return tOutputList;
 408
 409 }