PWGCF/FEMTOSCOPY/AliFemtoUser/AliFemtoBPLCMS3DCorrFctnEMCIC.cxx

   1 /***************************************************************************
   2  *
   3  * $Id: AliFemtoBPLCMS3DCorrFctnEMCIC.cxx  $
   4  *
   5  * Author: Nicolas Bock, Ohio State University, bock@mps.ohio-state.edu
   6  ***************************************************************************
   7  *
   8  * Description: Calculates of the 3D Correlation Function, and also
   9  *              produces histograms to calculate Energy Momentum Conservation
  10  *              Induced Correlations  (EMCICs)
  11  *
  12  * This Class produces the following histograms as function of Qinv
  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  * The class is derived from AliFemtoBPLCMS3DCorrFctn, therefore it produces
  20  * also the histograms in that class.
  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 //////////////////////////////////////////////////////////////////////////
  31 //                                                                       //
  32 // AliFemtoBPLCMS3DCorrFctn: a class to calculate 3D correlation         //
  33 // for pairs of identical particles.                                     //
  34 // It also stored the weighted qinv per bin histogram for the coulomb    //
  35 // correction.                                                           //
  36 // In analysis the function should be first created in a macro, then     //
  37 // added to the analysis, and at the end of the macro the procedure to   //
  38 // write out histograms should be called.                                //
  39 //                                                                       //
  40 ///////////////////////////////////////////////////////////////////////////
  41
  42 #include "AliFemtoBPLCMS3DCorrFctnEMCIC.h"
  43 #include "AliFemtoKTPairCut.h"
  44 #include "AliFemtoAnalysisReactionPlane.h"
  45 //#include "AliFemtoHisto.h"
  46 #include <cstdio>
  47 #include <TVector2.h>
  48
  49 #ifdef __ROOT__
  50 ClassImp(AliFemtoBPLCMS3DCorrFctnEMCIC)
  51 #endif
  52
  53 //____________________________
  54 AliFemtoBPLCMS3DCorrFctnEMCIC::AliFemtoBPLCMS3DCorrFctnEMCIC(char* title, const int& nbins, const float& QLo, const float& QHi)
  55 :
  56 AliFemtoCorrFctn(),
  57 //fEnergyTotalReal(0),
  58 // fEnergyMultReal(0),
  59 //fPzMultReal(0),
  60 //fPtMultReal(0),
  61   fNumerator(0),
  62   fDenominator(0),
  63   fEnergyTotalMix(0),
  64   fEnergyMultMix(0),
  65   fPzMultMix(0),
  66   fPtMultMix(0),
  67   fUseRPSelection(0)
  68 {
  69
  70   // set up numerator
  71   char tTitNum[101] = "Num";
  72   strncat(tTitNum,title, 100);
  73   fNumerator = new TH3D(tTitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  74   // set up denominator
  75   char tTitDen[101] = "Den";
  76   strncat(tTitDen,title, 100);
  77   fDenominator = new TH3D(tTitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  78
  79   //Setup EnergyTotalReal
  80   /*char tTitNum1[100] = "ESumReal";
  81   strncat(tTitNum1,title, 100);
  82   fEnergyTotalReal = new TH3D(tTitNum1,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  83
  84   //Setup EnergyMultReal
  85   char tTitNum2[100] = "EMultReal";
  86   strncat(tTitNum2,title, 100);
  87   fEnergyMultReal = new TH3D(tTitNum2,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  88
  89   //Setup Pz MultReal
  90   char tTitNum3[100] = "PzMultReal";
  91   strncat(tTitNum3,title, 100);
  92   fPzMultReal = new TH3D(tTitNum3,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  93
  94   //Setup Pt MultReal
  95   char tTitNum4[100] = "PtMultReal";
  96   strncat(tTitNum4,title, 100);
  97   fPtMultReal = new TH3D(tTitNum4,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);  */
  98
  99   //Setup EnergyTotalMix
 100   char tTitNum5[101] = "ESumMix";
 101   strncat(tTitNum5,title, 100);
 102   fEnergyTotalMix = new TH3D(tTitNum5,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 103
 104   //Setup EnergyMultMix
 105   char tTitNum6[101] = "EMultMix";
 106   strncat(tTitNum6,title, 100);
 107   fEnergyMultMix = new TH3D(tTitNum6,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 108
 109   //Setup Pz MultMix
 110   char tTitNum7[101] = "PzMultMix";
 111   strncat(tTitNum7,title, 100);
 112   fPzMultMix = new TH3D(tTitNum7,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 113
 114   //Setup Pt MultMix
 115   char tTitNum8[101] = "PtMultMix";
 116   strncat(tTitNum8,title, 100);
 117   fPtMultMix = new TH3D(tTitNum8,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
 118   // To enable error bar calculation
 119
 120   /*fEnergyTotalReal->Sumw2();
 121   fEnergyMultReal->Sumw2();
 122   fPzMultReal->Sumw2();
 123   fPtMultReal->Sumw2();  */
 124   fNumerator->Sumw2();
 125   fDenominator->Sumw2();
 126   fEnergyTotalMix->Sumw2();
 127   fEnergyMultMix->Sumw2();
 128   fPzMultMix->Sumw2();
 129   fPtMultMix->Sumw2();
 130
 131 }
 132
 133
 134 // Variable bin size constructor :
 135 //qBins array of low-edges for each bin. This is an array of size nbins+1
 136 AliFemtoBPLCMS3DCorrFctnEMCIC::AliFemtoBPLCMS3DCorrFctnEMCIC(char* title, const int& nbins, const float* qBins):
 137 AliFemtoCorrFctn(),
 138   fNumerator(0),
 139   fDenominator(0),
 140   fEnergyTotalMix(0),
 141   fEnergyMultMix(0),
 142   fPzMultMix(0),
 143   fPtMultMix(0),
 144   fUseRPSelection(0)
 145 {
 146
 147   // set up numerator
 148   char tTitNum[101] = "Num";
 149   strncat(tTitNum,title, 100);
 150   fNumerator = new TH3D(tTitNum,title,nbins,qBins,nbins,qBins,nbins,qBins);
 151   // set up denominator
 152   char tTitDen[101] = "Den";
 153   strncat(tTitDen,title, 100);
 154   fDenominator = new TH3D(tTitDen,title,nbins,qBins,nbins,qBins,nbins,qBins);
 155
 156   //Setup EnergyTotalReal
 157   /*char tTitNum1[100] = "ESumReal";
 158   strncat(tTitNum1,title, 100);
 159   fEnergyTotalReal = new TH3D(tTitNum1,title,nbins,qBins,nbins,qBins,nbins,qBins);
 160
 161   //Setup EnergyMultReal
 162   char tTitNum2[100] = "EMultReal";
 163   strncat(tTitNum2,title, 100);
 164   fEnergyMultReal = new TH3D(tTitNum2,title,nbins,qBins,nbins,qBins,nbins,qBins);
 165
 166   //Setup Pz MultReal
 167   char tTitNum3[100] = "PzMultReal";
 168   strncat(tTitNum3,title, 100);
 169   fPzMultReal = new TH3D(tTitNum3,title,nbins,qBins,nbins,qBins,nbins,qBins);
 170
 171   //Setup Pt MultReal
 172   char tTitNum4[100] = "PtMultReal";
 173   strncat(tTitNum4,title, 100);
 174   fPtMultReal = new TH3D(tTitNum4,title,nbins,qBins,nbins,qBins,nbins,qBins);  */
 175
 176   //Setup EnergyTotalMix
 177   char tTitNum5[101] = "ESumMix";
 178   strncat(tTitNum5,title, 100);
 179   fEnergyTotalMix = new TH3D(tTitNum5,title,nbins,qBins,nbins,qBins,nbins,qBins);
 180
 181   //Setup EnergyMultMix
 182   char tTitNum6[101] = "EMultMix";
 183   strncat(tTitNum6,title, 100);
 184   fEnergyMultMix = new TH3D(tTitNum6,title,nbins,qBins,nbins,qBins,nbins,qBins);
 185
 186   //Setup Pz MultMix
 187   char tTitNum7[101] = "PzMultMix";
 188   strncat(tTitNum7,title, 100);
 189   fPzMultMix = new TH3D(tTitNum7,title,nbins,qBins,nbins,qBins,nbins,qBins);
 190
 191   //Setup Pt MultMix
 192   char tTitNum8[101] = "PtMultMix";
 193   strncat(tTitNum8,title, 100);
 194   fPtMultMix = new TH3D(tTitNum8,title,nbins,qBins,nbins,qBins,nbins,qBins);
 195   // To enable error bar calculation
 196
 197   /*fEnergyTotalReal->Sumw2();
 198   fEnergyMultReal->Sumw2();
 199   fPzMultReal->Sumw2();
 200   fPtMultReal->Sumw2();  */
 201   fNumerator->Sumw2();
 202   fDenominator->Sumw2();
 203   fEnergyTotalMix->Sumw2();
 204   fEnergyMultMix->Sumw2();
 205   fPzMultMix->Sumw2();
 206   fPtMultMix->Sumw2();
 207 }
 208
 209
 210
 211
 212
 213
 214
 215 AliFemtoBPLCMS3DCorrFctnEMCIC::AliFemtoBPLCMS3DCorrFctnEMCIC(const AliFemtoBPLCMS3DCorrFctnEMCIC& aCorrFctn) :
 216   AliFemtoCorrFctn(aCorrFctn),
 217   /*fEnergyTotalReal(0),
 218   fEnergyMultReal(0),
 219   fPzMultReal(0),
 220   fPtMultReal(0),*/
 221   fNumerator(0),
 222   fDenominator(0),
 223   fEnergyTotalMix (0),
 224   fEnergyMultMix (0),
 225   fPzMultMix(0),
 226   fPtMultMix(0),
 227   fUseRPSelection(0)
 228 {
 229   // Copy constructor
 230   fNumerator = new TH3D(*aCorrFctn.fNumerator);
 231   fDenominator = new TH3D(*aCorrFctn.fDenominator);
 232   /*  fEnergyTotalReal = new TH3D(*aCorrFctn.fEnergyTotalReal);
 233   fEnergyMultReal = new TH3D(*aCorrFctn.fEnergyMultReal);
 234   fPzMultReal = new TH3D(*aCorrFctn.fPzMultReal);
 235   fPtMultReal = new TH3D(*aCorrFctn.fPtMultReal);  */
 236   fEnergyTotalMix = new TH3D(*aCorrFctn.fEnergyTotalMix);
 237   fEnergyMultMix = new TH3D(*aCorrFctn.fEnergyMultMix);
 238   fPzMultMix = new TH3D(*aCorrFctn.fPzMultMix);
 239   fPtMultMix = new TH3D(*aCorrFctn.fPtMultMix);
 240 }
 241 //____________________________
 242 AliFemtoBPLCMS3DCorrFctnEMCIC::~AliFemtoBPLCMS3DCorrFctnEMCIC(){
 243   // Destructor
 244   /*  delete fEnergyTotalReal;
 245   delete fEnergyMultReal;
 246   delete fPzMultReal;
 247   delete fPtMultReal;  */
 248   delete fNumerator;
 249   delete fDenominator;
 250   delete fEnergyTotalMix;
 251   delete fEnergyMultMix;
 252   delete fPzMultMix;
 253   delete fPtMultMix;
 254 }
 255 //_________________________
 256 AliFemtoBPLCMS3DCorrFctnEMCIC& AliFemtoBPLCMS3DCorrFctnEMCIC::operator=(const AliFemtoBPLCMS3DCorrFctnEMCIC& aCorrFctn)
 257 {
 258   // assignment operator
 259   if (this == &aCorrFctn)
 260     return *this;
 261   if (fNumerator) delete fNumerator;
 262   fNumerator = new TH3D(*aCorrFctn.fNumerator);
 263   if (fDenominator) delete fDenominator;
 264   fDenominator = new TH3D(*aCorrFctn.fDenominator);
 265   //Emcics
 266   /*  if (fEnergyTotalReal) delete fEnergyTotalReal;
 267   fEnergyTotalReal = new TH3D(*aCorrFctn.fEnergyTotalReal);
 268   if (fEnergyMultReal) delete fEnergyMultReal;
 269   fEnergyMultReal = new TH3D(*aCorrFctn.fEnergyMultReal);
 270   if (fPzMultReal) delete fPzMultReal;
 271   fPzMultReal = new TH3D(*aCorrFctn.fPzMultReal);
 272   if (fPtMultReal) delete fPtMultReal;
 273   fPtMultReal = new TH3D(*aCorrFctn.fPtMultReal);  */
 274   if (fEnergyTotalMix) delete fEnergyTotalMix;
 275   fEnergyTotalMix = new TH3D(*aCorrFctn.fEnergyTotalMix);
 276   if (fEnergyMultMix) delete fEnergyMultMix;
 277   fEnergyMultMix = new TH3D(*aCorrFctn.fEnergyMultMix);
 278   if (fPzMultMix) delete fPzMultMix;
 279   fPzMultMix = new TH3D(*aCorrFctn.fPzMultMix);
 280   if (fPtMultMix) delete fPtMultMix;
 281   fPtMultMix = new TH3D(*aCorrFctn.fPtMultMix);
 282
 283   fUseRPSelection = aCorrFctn.fUseRPSelection;
 284
 285   return *this;
 286 }
 287
 288 //_________________________
 289 void AliFemtoBPLCMS3DCorrFctnEMCIC::WriteOutHistos(){
 290
 291   fNumerator->Write();
 292   fDenominator->Write();
 293   //fEnergyTotalReal->Write();
 294   //fEnergyMultReal->Write();
 295   //fPzMultReal->Write();
 296   //fPtMultReal->Write();
 297   fEnergyTotalMix->Write();
 298   fEnergyMultMix->Write();
 299   fPzMultMix->Write();
 300   fPtMultMix->Write();
 301   //cout << "write histos emcics" << endl;
 302 }
 303 //______________________________
 304 TList* AliFemtoBPLCMS3DCorrFctnEMCIC::GetOutputList()
 305 {
 306   // Prepare the list of objects to be written to the output
 307   TList *tOutputList = new TList();
 308   cout << "Getting list from CFemcic" << endl;
 309   tOutputList->Add(fNumerator);
 310   tOutputList->Add(fDenominator);
 311   /*tOutputList->Add(fEnergyTotalReal);
 312   tOutputList->Add(fEnergyMultReal);
 313   tOutputList->Add(fPzMultReal);
 314   tOutputList->Add(fPtMultReal);     */
 315   tOutputList->Add(fEnergyTotalMix );
 316   tOutputList->Add(fEnergyMultMix );
 317   tOutputList->Add(fPzMultMix);
 318   tOutputList->Add(fPtMultMix);
 319   return tOutputList;
 320 }
 321
 322
 323
 324 //____________________________
 325 void AliFemtoBPLCMS3DCorrFctnEMCIC::AddRealPair( AliFemtoPair* pair){
 326   // perform operations on real pairs
 327
 328   if (fPairCut){
 329     if (fUseRPSelection) {
 330       AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
 331       if (!ktc) {
 332         cout << "RP aware cut requested, but not connected to the CF" << endl;
 333         if (!(fPairCut->Pass(pair))) return;
 334       }
 335       else {
 336         AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
 337         if (!arp) {
 338           cout << "RP aware cut requested, but not connected to the CF" << endl;
 339           if (!(fPairCut->Pass(pair))) return;
 340         }
 341         else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
 342       }
 343     }
 344     else
 345       if (!(fPairCut->Pass(pair))) return;
 346   }
 347
 348   double qOut  = fabs(pair->QOutCMS());
 349   double qSide = fabs( pair->QSideCMS());
 350   double qLong = fabs(pair->QLongCMS());
 351
 352   fNumerator->Fill(qOut,qSide,qLong);
 353
 354   /*AliFemtoLorentzVector tMom1 = pair->Track1()->FourMomentum();
 355   AliFemtoLorentzVector tMom2 = pair->Track2()->FourMomentum();
 356   double tE1 = tMom1.e();
 357   double tE2 = tMom2.e();
 358   double tPz1 = tMom1.pz();
 359   double tPz2 = tMom2.pz();
 360   TVector2 tPt1;
 361   TVector2 tPt2;
 362   tPt1.Set(tMom1.px(),tMom1.py());
 363   tPt2.Set(tMom2.px(),tMom2.py());
 364
 365
 366   double tPt1DotPt2 = tPt1*tPt2;
 367
 368   fEnergyTotalReal->Fill(qOut,qSide,qLong,tE1+tE2);
 369   fEnergyMultReal->Fill(qOut,qSide,qLong,tE1*tE2);
 370   fPzMultReal->Fill(qOut,qSide,qLong,tPz1*tPz2);
 371   fPtMultReal->Fill(qOut,qSide,qLong,tPt1DotPt2); */
 372 }
 373
 374
 375 //____________________________
 376 void AliFemtoBPLCMS3DCorrFctnEMCIC::AddMixedPair( AliFemtoPair* pair){
 377   // perform operations on mixed pairs
 378 //   if (fPairCut){
 379 //     if (!(fPairCut->Pass(pair))) return;
 380 //   }
 381   if (fPairCut){
 382     if (fUseRPSelection) {
 383       AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
 384       if (!ktc) {
 385         cout << "RP aware cut requested, but not connected to the CF" << endl;
 386         if (!(fPairCut->Pass(pair))) return;
 387       }
 388       else {
 389         AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
 390         if (!arp) {
 391           cout << "RP aware cut requested, but not connected to the CF" << endl;
 392           if (!(fPairCut->Pass(pair))) return;
 393         }
 394         else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
 395       }
 396     }
 397     else
 398       if (!(fPairCut->Pass(pair))) return;
 399   }
 400
 401
 402   double qOut = fabs(pair->QOutCMS());
 403   double qSide = fabs(pair->QSideCMS());
 404   double qLong = fabs(pair->QLongCMS());
 405
 406   fDenominator->Fill(qOut,qSide,qLong);
 407
 408   AliFemtoLorentzVector tMom1 = pair->Track1()->FourMomentum();
 409   AliFemtoLorentzVector tMom2 = pair->Track2()->FourMomentum();
 410   double tE1 = tMom1.e();
 411   double tE2 = tMom2.e();
 412   double tPz1 = tMom1.pz();
 413   double tPz2 = tMom2.pz();
 414   TVector2 tPt1;
 415   TVector2 tPt2;
 416   tPt1.Set(tMom1.px(),tMom1.py());
 417   tPt2.Set(tMom2.px(),tMom2.py());
 418   double tPt1DotPt2 = tPt1*tPt2;
 419
 420   fEnergyTotalMix->Fill(qOut,qSide,qLong,tE1+tE2);
 421   fEnergyMultMix->Fill(qOut,qSide,qLong,tE1*tE2);
 422   fPzMultMix->Fill(qOut,qSide,qLong,tPz1*tPz2);
 423   fPtMultMix->Fill(qOut,qSide,qLong,tPt1DotPt2);
 424
 425 }
 426
 427
 428 void AliFemtoBPLCMS3DCorrFctnEMCIC::SetUseRPSelection(unsigned short aRPSel)
 429 {
 430   fUseRPSelection = aRPSel;
 431 }