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