PWGCF/FEMTOSCOPY/AliFemto/AliFemtoQinvCorrFctn.cxx

   1 ///////////////////////////////////////////////////////////////////////////
   2 //                                                                       //
   3 // AliFemtoQinvCorrFctn:                                                 //
   4 // a simple Q-invariant correlation function                             //
   5 //                                                                       //
   6 ///////////////////////////////////////////////////////////////////////////
   7
   8 #include "AliFemtoQinvCorrFctn.h"
   9 //#include "AliFemtoHisto.h"
  10 #include <cstdio>
  11
  12 #ifdef __ROOT__
  13 ClassImp(AliFemtoQinvCorrFctn)
  14 #endif
  15
  16 //____________________________
  17 AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
  18   fNumerator(0),
  19   fDenominator(0),
  20   fRatio(0),
  21   fkTMonitor(0),
  22   fDetaDphiscal(kFALSE),
  23   fPairKinematics(kFALSE),
  24   fRaddedps(1.2),
  25   fNumDEtaDPhiS(0),
  26   fDenDEtaDPhiS(0),
  27   PairReader(0)// ,
  28   // fTrack1(NULL),
  29   // fTrack2(NULL)
  30
  31 {
  32   // set up numerator
  33   //  title = "Num Qinv (MeV/c)";
  34   char tTitNum[101] = "Num";
  35   strncat(tTitNum,title, 100);
  36   fNumerator = new TH1D(tTitNum,title,nbins,QinvLo,QinvHi);
  37   // set up denominator
  38   //title = "Den Qinv (MeV/c)";
  39   char tTitDen[101] = "Den";
  40   strncat(tTitDen,title, 100);
  41   fDenominator = new TH1D(tTitDen,title,nbins,QinvLo,QinvHi);
  42   // set up ratio
  43   //title = "Ratio Qinv (MeV/c)";
  44   char tTitRat[101] = "Rat";
  45   strncat(tTitRat,title, 100);
  46   fRatio = new TH1D(tTitRat,title,nbins,QinvLo,QinvHi);
  47
  48   char tTitkT[101] = "kTDep";
  49   strncat(tTitkT,title, 100);
  50   fkTMonitor = new TH1D(tTitkT,title,250,0.0,5.0);
  51
  52   char tTitNumDeDp[101] = "NumDEtaDPhiS";
  53   strncat(tTitNumDeDp,title, 100);
  54   fNumDEtaDPhiS = new TH2D(tTitNumDeDp,title,500,-0.2*TMath::Pi(),0.2*TMath::Pi(),500,-0.5,0.5);
  55
  56   char tTitDenDeDp[101] = "DenDEtaDPhiS";
  57   strncat(tTitDenDeDp,title, 100);
  58   fDenDEtaDPhiS = new TH2D(tTitDenDeDp,title,500,-0.2*TMath::Pi(),0.2*TMath::Pi(),500,-0.5,0.5);
  59
  60   char tTitPair[101] = "Pair";
  61   strncat(tTitPair,title, 100);
  62   PairReader = new TNtuple(tTitPair,title,  "px1:py1:pz1:e1:px2:py2:pz2:e2");
  63
  64
  65   // this next bit is unfortunately needed so that we can have many histos of same "title"
  66   // it is neccessary if we typedef TH1D to TH1d (which we do)
  67   //fNumerator->SetDirectory(0);
  68   //fDenominator->SetDirectory(0);
  69   //fRatio->SetDirectory(0);
  70
  71   // to enable error bar calculation...
  72   fNumerator->Sumw2();
  73   fDenominator->Sumw2();
  74   fRatio->Sumw2();
  75   fkTMonitor->Sumw2();
  76
  77   fNumDEtaDPhiS->Sumw2();
  78   fDenDEtaDPhiS->Sumw2();
  79
  80 }
  81
  82 //____________________________
  83 AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(const AliFemtoQinvCorrFctn& aCorrFctn) :
  84   AliFemtoCorrFctn(),
  85   fNumerator(0),
  86   fDenominator(0),
  87   fRatio(0),
  88   fkTMonitor(0),
  89   fDetaDphiscal(kFALSE),
  90   fPairKinematics(kFALSE),
  91   fRaddedps(1.2),
  92   fNumDEtaDPhiS(0),
  93   fDenDEtaDPhiS(0),
  94   PairReader(0)// ,
  95   // fTrack1(NULL),
  96   // fTrack2(NULL)
  97
  98 {
  99   // copy constructor
 100   fNumerator = new TH1D(*aCorrFctn.fNumerator);
 101   fDenominator = new TH1D(*aCorrFctn.fDenominator);
 102   fRatio = new TH1D(*aCorrFctn.fRatio);
 103   fkTMonitor = new TH1D(*aCorrFctn.fkTMonitor);
 104
 105   fNumDEtaDPhiS = new TH2D(*aCorrFctn.fNumDEtaDPhiS);
 106   fDenDEtaDPhiS = new TH2D(*aCorrFctn.fDenDEtaDPhiS);
 107
 108   fDetaDphiscal = aCorrFctn.fDetaDphiscal;
 109   fRaddedps = aCorrFctn.fRaddedps;
 110
 111   fPairKinematics = aCorrFctn.fPairKinematics;
 112
 113   if (aCorrFctn.PairReader)
 114     PairReader = (TNtuple*)aCorrFctn.PairReader;
 115
 116 }
 117 //____________________________
 118 AliFemtoQinvCorrFctn::~AliFemtoQinvCorrFctn(){
 119   // destructor
 120   delete fNumerator;
 121   delete fDenominator;
 122   delete fRatio;
 123   delete fkTMonitor;
 124   delete fNumDEtaDPhiS;
 125   delete fDenDEtaDPhiS;
 126   delete PairReader;
 127
 128 }
 129 //_________________________
 130 AliFemtoQinvCorrFctn& AliFemtoQinvCorrFctn::operator=(const AliFemtoQinvCorrFctn& aCorrFctn)
 131 {
 132   // assignment operator
 133   if (this == &aCorrFctn)
 134     return *this;
 135
 136   if (fNumerator) delete fNumerator;
 137   fNumerator = new TH1D(*aCorrFctn.fNumerator);
 138   if (fDenominator) delete fDenominator;
 139   fDenominator = new TH1D(*aCorrFctn.fDenominator);
 140   if (fRatio) delete fRatio;
 141   fRatio = new TH1D(*aCorrFctn.fRatio);
 142   if (fkTMonitor) delete fkTMonitor;
 143   fkTMonitor = new TH1D(*aCorrFctn.fkTMonitor);
 144
 145   if (fNumDEtaDPhiS) delete fNumDEtaDPhiS;
 146   fNumDEtaDPhiS = new TH2D(*aCorrFctn.fNumDEtaDPhiS);
 147   if (fDenDEtaDPhiS) delete fDenDEtaDPhiS;
 148   fDenDEtaDPhiS = new TH2D(*aCorrFctn.fDenDEtaDPhiS);
 149
 150   fDetaDphiscal = aCorrFctn.fDetaDphiscal;
 151   fRaddedps = aCorrFctn.fRaddedps;
 152
 153   fPairKinematics = aCorrFctn.fPairKinematics;
 154
 155   if (aCorrFctn.PairReader)
 156     PairReader = (TNtuple*)aCorrFctn.PairReader;
 157
 158   return *this;
 159 }
 160
 161 //_________________________
 162 void AliFemtoQinvCorrFctn::Finish(){
 163   // here is where we should normalize, fit, etc...
 164   // we should NOT Draw() the histos (as I had done it below),
 165   // since we want to insulate ourselves from root at this level
 166   // of the code.  Do it instead at root command line with browser.
 167   //  fNumerator->Draw();
 168   //fDenominator->Draw();
 169   //fRatio->Draw();
 170   fRatio->Divide(fNumerator,fDenominator,1.0,1.0);
 171
 172 }
 173
 174 //____________________________
 175 AliFemtoString AliFemtoQinvCorrFctn::Report(){
 176   // construct report
 177   string stemp = "Qinv Correlation Function Report:\n";
 178   char ctemp[100];
 179   snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
 180   stemp += ctemp;
 181   snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
 182   stemp += ctemp;
 183   snprintf(ctemp , 100, "Number of entries in ratio:\t%E\n",fRatio->GetEntries());
 184   stemp += ctemp;
 185   //  stemp += mCoulombWeight->Report();
 186   AliFemtoString returnThis = stemp;
 187   return returnThis;
 188 }
 189 //____________________________
 190 void AliFemtoQinvCorrFctn::AddRealPair(AliFemtoPair* pair){
 191   // add true pair
 192   if (fPairCut)
 193     if (!fPairCut->Pass(pair)) return;
 194
 195   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 196
 197   fNumerator->Fill(tQinv);
 198   fkTMonitor->Fill(pair->KT());
 199
 200
 201 //_______________________________________
 202   if (fDetaDphiscal) {
 203
 204     double phi1 = pair->Track1()->Track()->P().Phi();
 205     double phi2 = pair->Track2()->Track()->P().Phi();
 206     double chg1 = pair->Track1()->Track()->Charge();
 207     double chg2 = pair->Track2()->Track()->Charge();
 208     double ptv1 = pair->Track1()->Track()->Pt();
 209     double ptv2 = pair->Track2()->Track()->Pt();
 210     double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
 211     double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
 212
 213     AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 214     Int_t magsign = 0;
 215
 216     if (!aodH) {
 217       //AliWarning("Could not get AODInputHandler");
 218     }
 219     else {
 220       AliAODEvent *fAOD // = new AliAODEvent()
 221         ;
 222       magsign = fAOD->GetMagneticField();
 223       fAOD = aodH->GetEvent();
 224     }
 225
 226
 227     Int_t fMagSign;
 228
 229     if (magsign > 1)
 230       fMagSign = 1;
 231     else if ( magsign < 1)
 232       fMagSign = -1;
 233     else
 234       fMagSign = magsign;
 235
 236     Double_t rad = fRaddedps;
 237
 238     double afsi0b = 0.07510020733*chg1*fMagSign*rad/ptv1;
 239     double afsi1b = 0.07510020733*chg2*fMagSign*rad/ptv2;
 240     Double_t dps6 =  phi2 - phi1 + TMath::ASin(afsi1b) - TMath::ASin(afsi0b);
 241     dps6 = TVector2::Phi_mpi_pi(dps6);
 242
 243     // Double_t dps = (phi1-phi2+(TMath::ASin(-0.075*chg1*fMagSign*rad/ptv1))-(TMath::ASin(-0.075*chg2*fMagSign*rad/ptv2)));
 244     // dps = TVector2::Phi_mpi_pi(dps);
 245     double etad = eta2 - eta1;
 246
 247     fNumDEtaDPhiS->Fill(dps6,etad);
 248   }
 249 //_______________________________________________________________
 250
 251   //  cout << "AliFemtoQinvCorrFctn::AddRealPair : " << pair->qInv() << " " << tQinv <<
 252   //" " << pair->track1().FourMomentum() << " " << pair->track2().FourMomentum() << endl;
 253 }
 254
 255 //____________________________
 256 void AliFemtoQinvCorrFctn::AddMixedPair(AliFemtoPair* pair){
 257   // add mixed (background) pair
 258   if (fPairCut)
 259     if (!fPairCut->Pass(pair)) return;
 260
 261   double weight = 1.0;
 262   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 263   fDenominator->Fill(tQinv,weight);
 264
 265   if (fPairKinematics) {
 266     AliFemtoParticle* fTrack1 = pair->Track1();
 267     AliFemtoParticle* fTrack2 = pair->Track2();
 268
 269     double px1 = fTrack1->FourMomentum().vect().x();
 270     double py1 = fTrack1->FourMomentum().vect().y();
 271     double pz1 = fTrack1->FourMomentum().vect().z();
 272     double e1 = fTrack1->FourMomentum().e();
 273
 274     double px2 = fTrack2->FourMomentum().vect().x();
 275     double py2 = fTrack2->FourMomentum().vect().y();
 276     double pz2 = fTrack2->FourMomentum().vect().z();
 277     double e2 = fTrack2->FourMomentum().e();
 278     PairReader->Fill(px1, py1, pz1, e1, px2, py2, pz2, e2);
 279   }
 280
 281 //_______________________________________
 282   if (fDetaDphiscal) {
 283
 284     double phi1 = pair->Track1()->Track()->P().Phi();
 285     double phi2 = pair->Track2()->Track()->P().Phi();
 286     double chg1 = pair->Track1()->Track()->Charge();
 287     double chg2 = pair->Track2()->Track()->Charge();
 288     double ptv1 = pair->Track1()->Track()->Pt();
 289     double ptv2 = pair->Track2()->Track()->Pt();
 290     double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
 291     double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
 292
 293     AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 294     Int_t magsign = 0;
 295
 296
 297     if (!aodH) {
 298       //AliWarning("Could not get AODInputHandler");
 299       cout << "Could not get AODInputHandler" << endl;
 300     }
 301     else {
 302       AliAODEvent *fAOD;
 303       fAOD = aodH->GetEvent();
 304       magsign = fAOD->GetMagneticField();
 305     }
 306
 307
 308     Int_t fMagSign;
 309     if (magsign > 1)
 310       fMagSign = 1;
 311     else if ( magsign < 1)
 312       fMagSign = -1;
 313     else
 314       fMagSign = magsign;
 315
 316     Double_t rad = fRaddedps;
 317
 318     double afsi0b = 0.07510020733*chg1*fMagSign*rad/ptv1;
 319     double afsi1b = 0.07510020733*chg2*fMagSign*rad/ptv2;
 320     Double_t dps6 =  phi2 - phi1 + TMath::ASin(afsi1b) - TMath::ASin(afsi0b);
 321     dps6 = TVector2::Phi_mpi_pi(dps6);
 322     double etad = eta2 - eta1;
 323
 324     fDenDEtaDPhiS->Fill(dps6,etad);
 325   }
 326 //_______________________________________________________________
 327
 328 }
 329 //____________________________
 330 void AliFemtoQinvCorrFctn::Write(){
 331   // Write out neccessary objects
 332   fNumerator->Write();
 333   fDenominator->Write();
 334   fkTMonitor->Write();
 335   if (fDetaDphiscal) {
 336     fNumDEtaDPhiS->Write();
 337     fDenDEtaDPhiS->Write();
 338   }
 339   if (fPairKinematics) {
 340     PairReader->Write();
 341   }
 342 }
 343 //______________________________
 344 TList* AliFemtoQinvCorrFctn::GetOutputList()
 345 {
 346   // Prepare the list of objects to be written to the output
 347   TList *tOutputList = new TList();
 348
 349   tOutputList->Add(fNumerator);
 350   tOutputList->Add(fDenominator);
 351   tOutputList->Add(fkTMonitor);
 352   if (fDetaDphiscal) {
 353     tOutputList->Add(fNumDEtaDPhiS);
 354     tOutputList->Add(fDenDEtaDPhiS);
 355   }
 356   if (fPairKinematics) {
 357     tOutputList->Add(PairReader);
 358   }
 359   return tOutputList;
 360 }
 361
 362 void AliFemtoQinvCorrFctn::CalculateDetaDphis(Bool_t dedpsc, Double_t rad) {
 363   fDetaDphiscal = dedpsc;
 364   fRaddedps = rad;
 365 }
 366
 367 void AliFemtoQinvCorrFctn::CalculatePairKinematics(Bool_t pk) {
 368   fPairKinematics = pk;
 369 }