PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctn3DSphericalEMCIC.cxx

   1 /***************************************************************************
   2  *
   3  * $Id: AliFemtoCorrFctn3DSphercicalEMCIC.h  $
   4  *
   5  * Author: Nicolas Bock, Ohio State University, bock@mps.ohio-state.edu
   6  ***************************************************************************
   7  *
   8  * Description: Calculates of the 3D Correlation Function in Spherical
   9  *              coordinates, and also produces histograms to calculate
  10  *              Energy Momentum Conservation Induced Correlations  (EMCICs)
  11  *
  12  * This Class produces the following histograms as function of Q, theta, phi
  13  * (for both real and mixed pairs):
  14  *        1)   E1 + E2
  15  *        2)   E1 * E2
  16  *        3)   Pt1*Pt2
  17  *        4)   Pz1*Pz2
  18  *
  19  * This class is similar to AliFemtoCorrFctn3DSpherical, but it uses Q
  20  * instead of K to do the binning.
  21  *
  22  * NOTE: The EMCIC histograms are not averaged in this class, to obtain
  23  * the average, the user needs to divide the real pair histograms by
  24  * the numerator, and the mixed pair histograms by the denominator.
  25  *
  26  ***************************************************************************
  27  *
  28  **************************************************************************/
  29
  30 #include "AliFemtoCorrFctn3DSphericalEMCIC.h"
  31 #include <TMath.h>
  32 #include <TVector2.h>
  33 #include <cstdio>
  34
  35 #ifdef __ROOT__
  36 ClassImp(AliFemtoCorrFctn3DSphericalEMCIC)
  37 #endif
  38
  39
  40 //____________________________
  41 AliFemtoCorrFctn3DSphericalEMCIC::AliFemtoCorrFctn3DSphericalEMCIC(char* title, const int& nqbins, const float& QLo, const float& QHi, const int& nphibins, const int& ncthetabins):
  42   AliFemtoCorrFctn(),
  43   fNumerator(0),
  44   fDenominator(0),
  45   fEnergyTotalReal(0),
  46   fEnergyMultReal(0),
  47   fPzMultReal(0),
  48   fPtMultReal(0),
  49   fEnergyTotalMix (0),
  50   fEnergyMultMix (0),
  51   fPzMultMix(0),
  52   fPtMultMix(0),
  53   fPairCut(0x0)
  54 {
  55
  56   // To better sample on phi shift bin low edge by binsize/2 = Pi/numBins.
  57   Double_t shiftPhi=TMath::Pi()/nphibins;
  58   // set up numerator
  59   char tTitNum[100] = "Num";
  60   strcat(tTitNum,title);
  61   fNumerator = new TH3D(tTitNum,title,nqbins,QLo,QHi,nphibins,
  62                         -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
  63   // set up denominator
  64   char tTitDen[100] = "Den";
  65   strcat(tTitDen,title);
  66   fDenominator = new TH3D(tTitDen,title,nqbins,QLo,QHi,nphibins,
  67                           -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
  68
  69   //Added histograms to calculate EMCICs , Nicolas Bock 19.01.2010
  70   //Setup EnergyTotalReal
  71   char tTitNum1[100] = "ESumReal";
  72   strcat(tTitNum1,title);
  73   fEnergyTotalReal = new TH3D(tTitNum1,title,nqbins,QLo,QHi,nphibins,
  74                               -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
  75
  76   //Setup EnergyMultReal
  77   char tTitNum2[100] = "EMultReal";
  78   strcat(tTitNum2,title);
  79   fEnergyMultReal = new TH3D(tTitNum2,title,nqbins,QLo,QHi,nphibins,
  80                              -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
  81
  82   //Setup Pz MultReal
  83   char tTitNum3[100] = "PzMultReal";
  84   strcat(tTitNum3,title);
  85   fPzMultReal = new TH3D(tTitNum3,title,nqbins,QLo,QHi,nphibins,
  86                          -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
  87
  88   //Setup Pt MultReal
  89   char tTitNum4[100] = "PtMultReal";
  90   strcat(tTitNum4,title);
  91   fPtMultReal = new TH3D(tTitNum4,title,nqbins,QLo,QHi,nphibins,
  92                          -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
  93
  94
  95
  96   //Setup EnergyTotalMix
  97   char tTitNum5[100] = "ESumMix";
  98   strcat(tTitNum5,title);
  99   fEnergyTotalMix = new TH3D(tTitNum5,title,nqbins,QLo,QHi,nphibins,
 100                              -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
 101
 102   //Setup EnergyMultMix
 103   char tTitNum6[100] = "EMultMix";
 104   strcat(tTitNum6,title);
 105   fEnergyMultMix = new TH3D(tTitNum6,title,nqbins,QLo,QHi,nphibins,
 106                             -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
 107
 108   //Setup Pz MultMix
 109   char tTitNum7[100] = "PzMultMix";
 110   strcat(tTitNum7,title);
 111   fPzMultMix = new TH3D(tTitNum7,title,nqbins,QLo,QHi,nphibins,
 112                         -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
 113
 114   //Setup Pt MultMix
 115   char tTitNum8[100] = "PtMultMix";
 116   strcat(tTitNum8,title);
 117   fPtMultMix = new TH3D(tTitNum8,title,nqbins,QLo,QHi,nphibins,
 118                         -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,
 119                         ncthetabins,-1.0,1.0);
 120
 121   // to enable error bar calculation...
 122   fNumerator->Sumw2();
 123   fDenominator->Sumw2();
 124   fEnergyTotalReal->Sumw2();
 125   fEnergyMultReal->Sumw2();
 126   fPzMultReal->Sumw2();
 127   fPtMultReal->Sumw2();
 128   fEnergyTotalMix->Sumw2();
 129   fEnergyMultMix->Sumw2();
 130   fPzMultMix->Sumw2();
 131   fPtMultMix->Sumw2();
 132
 133
 134 }
 135
 136 AliFemtoCorrFctn3DSphericalEMCIC::AliFemtoCorrFctn3DSphericalEMCIC(const AliFemtoCorrFctn3DSphericalEMCIC& aCorrFctn) :
 137   AliFemtoCorrFctn(),
 138   fNumerator(0),
 139   fDenominator(0),
 140   fEnergyTotalReal(0),
 141   fEnergyMultReal(0),
 142   fPzMultReal(0),
 143   fPtMultReal(0),
 144   fEnergyTotalMix (0),
 145   fEnergyMultMix (0),
 146   fPzMultMix(0),
 147   fPtMultMix(0),
 148   fPairCut(0x0)
 149 {
 150   // Copy constructor
 151   fNumerator = new TH3D(*aCorrFctn.fNumerator);
 152   fDenominator = new TH3D(*aCorrFctn.fDenominator);
 153   fEnergyTotalReal = new TH3D(*aCorrFctn.fEnergyTotalReal);
 154   fEnergyMultReal = new TH3D(*aCorrFctn.fEnergyMultReal);
 155   fPzMultReal = new TH3D(*aCorrFctn.fPzMultReal);
 156   fPtMultReal = new TH3D(*aCorrFctn.fPtMultReal);
 157   fEnergyTotalMix = new TH3D(*aCorrFctn.fEnergyTotalMix);
 158   fEnergyMultMix = new TH3D(*aCorrFctn.fEnergyMultMix);
 159   fPzMultMix = new TH3D(*aCorrFctn.fPzMultMix);
 160   fPtMultMix = new TH3D(*aCorrFctn.fPtMultMix);
 161   fPairCut = aCorrFctn.fPairCut;
 162 }
 163 //____________________________
 164 AliFemtoCorrFctn3DSphericalEMCIC::~AliFemtoCorrFctn3DSphericalEMCIC(){
 165   // Destructor
 166   delete fNumerator;
 167   delete fDenominator;
 168   delete fEnergyTotalReal;
 169   delete fEnergyMultReal;
 170   delete fPzMultReal;
 171   delete fPtMultReal;
 172   delete fEnergyTotalMix;
 173   delete fEnergyMultMix;
 174   delete fPzMultMix;
 175   delete fPtMultMix;
 176 }
 177 //_________________________
 178 AliFemtoCorrFctn3DSphericalEMCIC& AliFemtoCorrFctn3DSphericalEMCIC::operator=(const AliFemtoCorrFctn3DSphericalEMCIC& aCorrFctn)
 179 {
 180   // assignment operator
 181   if (this == &aCorrFctn)
 182     return *this;
 183
 184   if (fNumerator) delete fNumerator;
 185   fNumerator = new TH3D(*aCorrFctn.fNumerator);
 186   if (fDenominator) delete fDenominator;
 187   fDenominator = new TH3D(*aCorrFctn.fDenominator);
 188   if (fEnergyTotalReal) delete fEnergyTotalReal;
 189   fEnergyTotalReal = new TH3D(*aCorrFctn.fEnergyTotalReal);
 190   if (fEnergyMultReal) delete fEnergyMultReal;
 191   fEnergyMultReal = new TH3D(*aCorrFctn.fEnergyMultReal);
 192   if (fPzMultReal) delete fPzMultReal;
 193   fPzMultReal = new TH3D(*aCorrFctn.fPzMultReal);
 194   if (fPtMultReal) delete fPtMultReal;
 195   fPtMultReal = new TH3D(*aCorrFctn.fPtMultReal);
 196   if (fEnergyTotalMix) delete fEnergyTotalMix;
 197   fEnergyTotalMix = new TH3D(*aCorrFctn.fEnergyTotalMix);
 198   if (fEnergyMultMix) delete fEnergyMultMix;
 199   fEnergyMultMix = new TH3D(*aCorrFctn.fEnergyMultMix);
 200   if (fPzMultMix) delete fPzMultMix;
 201   fPzMultMix = new TH3D(*aCorrFctn.fPzMultMix);
 202   if (fPtMultMix) delete fPtMultMix;
 203   fPtMultMix = new TH3D(*aCorrFctn.fPtMultMix);
 204   fPairCut = aCorrFctn.fPairCut;
 205
 206   return *this;
 207 }
 208
 209 //_________________________
 210 void AliFemtoCorrFctn3DSphericalEMCIC::WriteOutHistos(){
 211   // Write out all histograms to file
 212   fNumerator->Write();
 213   fDenominator->Write();
 214   fEnergyTotalReal->Write();
 215   fEnergyMultReal->Write();
 216   fPzMultReal->Write();
 217   fPtMultReal->Write();
 218   fEnergyTotalMix->Write();
 219   fEnergyMultMix->Write();
 220   fPzMultMix->Write();
 221   fPtMultMix->Write();
 222 }
 223 //______________________________
 224 TList* AliFemtoCorrFctn3DSphericalEMCIC::GetOutputList()
 225 {
 226   // Prepare the list of objects to be written to the output
 227   TList *tOutputList = new TList();
 228
 229   tOutputList->Add(fNumerator);
 230   tOutputList->Add(fDenominator);
 231   tOutputList->Add(fEnergyTotalReal);
 232   tOutputList->Add(fEnergyMultReal);
 233   tOutputList->Add(fPzMultReal);
 234   tOutputList->Add(fPtMultReal);
 235   tOutputList->Add(fEnergyTotalMix );
 236   tOutputList->Add(fEnergyMultMix );
 237   tOutputList->Add(fPzMultMix);
 238   tOutputList->Add(fPtMultMix);
 239   return tOutputList;
 240 }
 241
 242 //_________________________
 243 void AliFemtoCorrFctn3DSphericalEMCIC::Finish(){
 244   // here is where we should normalize, fit, etc...
 245 }
 246
 247 //____________________________
 248 AliFemtoString AliFemtoCorrFctn3DSphericalEMCIC::Report(){
 249   // Construct the report
 250   string stemp = "PRF Frame SphericalEMCIC 3D Correlation Function Report:\n";
 251   char ctemp[100];
 252   sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
 253   stemp += ctemp;
 254   sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
 255   stemp += ctemp;
 256
 257   if (fPairCut){
 258     sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
 259     stemp += ctemp;
 260     stemp += fPairCut->Report();
 261   }
 262   else{
 263     sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
 264     stemp += ctemp;
 265   }
 266
 267   //
 268   AliFemtoString returnThis = stemp;
 269   return returnThis;
 270 }
 271 //____________________________
 272 void AliFemtoCorrFctn3DSphericalEMCIC::AddRealPair( AliFemtoPair* pair){
 273   // perform operations on real pairs
 274   if (fPairCut){
 275     if (!(fPairCut->Pass(pair))) return;
 276   }
 277
 278   //                          //Changed K to Q to be in LCMS, N. Bock
 279   double tKO = pair->QOutCMS();  //KOut();
 280   double tKS = pair->QSideCMS();  //KSide();
 281   double tKL = pair->QLongCMS();  //  KLong();
 282
 283   double tKR = sqrt(tKO*tKO + tKS*tKS + tKL*tKL);
 284   double tKC;
 285   if ( fabs(tKR) < 1e-10 ) tKC = 0.0;
 286   else tKC=tKL/tKR;
 287   double tKP=atan2(tKS,tKO);
 288
 289   fNumerator->Fill(tKR,tKP,tKC);
 290
 291   // EMCICs
 292   AliFemtoLorentzVector tMom1 = pair->Track1()->FourMomentum();
 293   AliFemtoLorentzVector tMom2 = pair->Track2()->FourMomentum();
 294   double tE1 = tMom1.e();
 295   double tE2 = tMom2.e();
 296   double tPz1 = tMom1.pz();
 297   double tPz2 = tMom2.pz();
 298   TVector2 tPt1;
 299   TVector2 tPt2;
 300   cout <<"AddRealPair" << "E1,E2, K,Ks,KL" <<
 301     tE1 <<" " << tE2<<" " <<  tKR <<" " <<  tKS <<" " << tKL<<endl;
 302   tPt1.Set(tMom1.px(),tMom1.py());
 303   tPt2.Set(tMom2.px(),tMom2.py());
 304   double tPt1DotPt2=tPt1*tPt2;
 305
 306   fEnergyTotalReal->Fill(tKR,tKP,tKC,tE1+tE2);
 307   fEnergyMultReal->Fill(tKR,tKP,tKC,tE1*tE2);
 308   fPzMultReal->Fill(tKR,tKP,tKC,tPz1*tPz2);
 309   fPtMultReal->Fill(tKR,tKP,tKC,tPt1DotPt2);
 310
 311
 312 }
 313 //____________________________
 314 void AliFemtoCorrFctn3DSphericalEMCIC::AddMixedPair( AliFemtoPair* pair){
 315   // perform operations on mixed pairs
 316   if (fPairCut){
 317     if (!(fPairCut->Pass(pair))) return;
 318   }
 319  //                          //Changed K to Q to be in LCMS, N. Bock
 320   double tKO = pair->QOutCMS();  //KOut();
 321   double tKS = pair->QSideCMS();  //KSide();
 322   double tKL = pair->QLongCMS();  //  KLong();
 323
 324   double tKR = sqrt(tKO*tKO + tKS*tKS + tKL*tKL);
 325   double tKC;
 326   if ( fabs(tKR) < 1e-10 ) tKC = 0.0;
 327   else tKC=tKL/tKR;
 328   double tKP=atan2(tKS,tKO);
 329
 330   fDenominator->Fill(tKR,tKP,tKC);
 331
 332   // EMCICs
 333   AliFemtoLorentzVector tMom1 = pair->Track1()->FourMomentum();
 334   AliFemtoLorentzVector tMom2 = pair->Track2()->FourMomentum();
 335   double tE1 = tMom1.e();
 336   double tE2 = tMom2.e();
 337   double tPz1 = tMom1.pz();
 338   double tPz2 = tMom2.pz();
 339   TVector2 tPt1;
 340   TVector2 tPt2;
 341   tPt1.Set(tMom1.px(),tMom1.py());
 342   tPt2.Set(tMom2.px(),tMom2.py());
 343   double tPt1DotPt2=tPt1*tPt2;
 344
 345   fEnergyTotalMix->Fill(tKR,tKP,tKC,tE1+tE2);
 346   fEnergyMultMix->Fill(tKR,tKP,tKC,tE1*tE2);
 347   fPzMultMix->Fill(tKR,tKP,tKC,tPz1*tPz2);
 348   fPtMultMix->Fill(tKR,tKP,tKC,tPt1DotPt2);
 349
 350 }