PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoTPCInnerCorrFctn.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 ///                                                                          ///
   3 /// AliFemtoTPCInnerCorrFctn - A correlation function that saves the         ///
   4 /// distance at the entrance to the TPC between two tracks as a function     ///
   5 /// of qinv                                                                  ///
   6 /// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
   7 ///                                                                          ///
   8 ////////////////////////////////////////////////////////////////////////////////
   9
  10 #include "AliFemtoTPCInnerCorrFctn.h"
  11 //#include "AliFemtoHisto.hh"
  12 #include <cstdio>
  13
  14 #ifdef __ROOT__
  15 ClassImp(AliFemtoTPCInnerCorrFctn)
  16 #endif
  17
  18 //____________________________
  19 AliFemtoTPCInnerCorrFctn::AliFemtoTPCInnerCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
  20   fDTPCNumerator(0),
  21   fDTPCDenominator(0),
  22   fRadDNumerator(0),
  23   fRadDDenominator(0),
  24   fRadius(100)
  25 {
  26   // set up numerator
  27   //  title = "Num Qinv (MeV/c)";
  28   char tTitNum[101] = "NumDTPC";
  29   strncat(tTitNum,title, 100);
  30   fDTPCNumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,100,0.0,20.0);
  31   // set up denominator
  32   //title = "Den Qinv (MeV/c)";
  33   char tTitDen[101] = "DenDTPC";
  34   strncat(tTitDen,title, 100);
  35   fDTPCDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,100,0.0,20.0);
  36
  37   char tTitNumR[101] = "NumRadD";
  38   strncat(tTitNumR,title, 100);
  39   fRadDNumerator = new TH2D(tTitNumR,title,50,-0.1,0.1,50,-0.1,0.1);
  40   // set up denominator
  41   //title = "Den Qinv (MeV/c)";
  42   char tTitDenR[101] = "DenRadD";
  43   strncat(tTitDenR,title, 100);
  44   fRadDDenominator = new TH2D(tTitDenR,title,50,-0.1,0.1,50,-0.1,0.1);
  45
  46   // to enable error bar calculation...
  47   fDTPCNumerator->Sumw2();
  48   fDTPCDenominator->Sumw2();
  49   fRadDNumerator->Sumw2();
  50   fRadDDenominator->Sumw2();
  51 }
  52
  53 //____________________________
  54 AliFemtoTPCInnerCorrFctn::AliFemtoTPCInnerCorrFctn(const AliFemtoTPCInnerCorrFctn& aCorrFctn) :
  55   AliFemtoCorrFctn(),
  56   fDTPCNumerator(0),
  57   fDTPCDenominator(0),
  58   fRadDNumerator(0),
  59   fRadDDenominator(0),
  60   fRadius(100)
  61 {
  62   // copy constructor
  63   if (aCorrFctn.fDTPCNumerator)
  64     fDTPCNumerator = new TH2D(*aCorrFctn.fDTPCNumerator);
  65   if (aCorrFctn.fDTPCDenominator)
  66     fDTPCDenominator = new TH2D(*aCorrFctn.fDTPCDenominator);
  67   if (aCorrFctn.fRadDNumerator)
  68     fRadDNumerator = new TH2D(*aCorrFctn.fRadDNumerator);
  69   if (aCorrFctn.fRadDDenominator)
  70     fRadDDenominator = new TH2D(*aCorrFctn.fRadDDenominator);
  71 }
  72 //____________________________
  73 AliFemtoTPCInnerCorrFctn::~AliFemtoTPCInnerCorrFctn(){
  74   // destructor
  75   delete fDTPCNumerator;
  76   delete fDTPCDenominator;
  77   delete fRadDNumerator;
  78   delete fRadDDenominator;
  79 }
  80 //_________________________
  81 AliFemtoTPCInnerCorrFctn& AliFemtoTPCInnerCorrFctn::operator=(const AliFemtoTPCInnerCorrFctn& aCorrFctn)
  82 {
  83   // assignment operator
  84   if (this == &aCorrFctn)
  85     return *this;
  86
  87   if (aCorrFctn.fDTPCNumerator)
  88     fDTPCNumerator = new TH2D(*aCorrFctn.fDTPCNumerator);
  89   else
  90     fDTPCNumerator = 0;
  91   if (aCorrFctn.fDTPCDenominator)
  92     fDTPCDenominator = new TH2D(*aCorrFctn.fDTPCDenominator);
  93   else
  94     fDTPCDenominator = 0;
  95
  96   if (aCorrFctn.fRadDNumerator)
  97     fRadDNumerator = new TH2D(*aCorrFctn.fRadDNumerator);
  98   else
  99     fRadDNumerator = 0;
 100   if (aCorrFctn.fRadDDenominator)
 101     fRadDDenominator = new TH2D(*aCorrFctn.fRadDDenominator);
 102   else
 103     fRadDDenominator = 0;
 104
 105   fRadius = aCorrFctn.fRadius;
 106
 107   return *this;
 108 }
 109 //_________________________
 110 void AliFemtoTPCInnerCorrFctn::Finish(){
 111   // here is where we should normalize, fit, etc...
 112   // we should NOT Draw() the histos (as I had done it below),
 113   // since we want to insulate ourselves from root at this level
 114   // of the code.  Do it instead at root command line with browser.
 115   //  mShareNumerator->Draw();
 116   //mShareDenominator->Draw();
 117   //mRatio->Draw();
 118
 119 }
 120
 121 //____________________________
 122 AliFemtoString AliFemtoTPCInnerCorrFctn::Report(){
 123   // create report
 124   string stemp = "Entrace TPC distance Correlation Function Report:\n";
 125   char ctemp[100];
 126   snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fDTPCNumerator->GetEntries());
 127   stemp += ctemp;
 128   snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDTPCDenominator->GetEntries());
 129   stemp += ctemp;
 130   //  stemp += mCoulombWeight->Report();
 131   AliFemtoString returnThis = stemp;
 132   return returnThis;
 133 }
 134 //____________________________
 135 void AliFemtoTPCInnerCorrFctn::AddRealPair( AliFemtoPair* pair){
 136   // add real (effect) pair
 137   if (fPairCut)
 138     if (!fPairCut->Pass(pair)) return;
 139
 140   double pih = TMath::Pi();
 141   double pit = TMath::Pi()*2;
 142
 143   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 144   double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
 145   double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
 146   double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
 147   double dist = sqrt(distx*distx + disty*disty + distz*distz);
 148
 149   fDTPCNumerator->Fill(tQinv, dist);
 150
 151   if (tQinv < 0.1) {
 152     double phi1 = pair->Track1()->Track()->P().Phi();
 153     double phi2 = pair->Track2()->Track()->P().Phi();
 154     double chg1 = pair->Track1()->Track()->Charge();
 155     double chg2 = pair->Track2()->Track()->Charge();
 156     double ptv1 = pair->Track1()->Track()->Pt();
 157     double ptv2 = pair->Track2()->Track()->Pt();
 158     double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
 159     double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
 160     double arg1 = -0.3 * 0.5 * chg1 * fRadius/(2*ptv1);
 161     double arg2 = -0.3 * 0.5 * chg2 * fRadius/(2*ptv2);
 162     double phid = phi2 - phi1 + TMath::ASin(arg2) - TMath::ASin(arg1);
 163
 164     while (phid>pih) phid -= pit;
 165     while (phid<-pih) phid += pit;
 166     //    dist = phi2 - phi1 + TMath::ASin(-0.3 * 0.5 * chg2 * fRadius/(2*ptv2)) - TMath::ASin(-0.3 * 0.5 * chg1 * fRadius/(2*ptv1));
 167     double etad = eta2 - eta1;
 168     fRadDNumerator->Fill(phid, etad);
 169   }
 170 }
 171 //____________________________
 172 void AliFemtoTPCInnerCorrFctn::AddMixedPair( AliFemtoPair* pair){
 173   // add mixed (background) pair
 174   if (fPairCut)
 175     if (!fPairCut->Pass(pair)) return;
 176
 177   double pih = TMath::Pi();
 178   double pit = TMath::Pi()*2;
 179
 180   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 181   double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
 182   double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
 183   double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
 184   double dist = sqrt(distx*distx + disty*disty + distz*distz);
 185
 186   fDTPCDenominator->Fill(tQinv,dist);
 187
 188   if (tQinv < 0.1) {
 189     double phi1 = pair->Track1()->Track()->P().Phi();
 190     double phi2 = pair->Track2()->Track()->P().Phi();
 191     double chg1 = pair->Track1()->Track()->Charge();
 192     double chg2 = pair->Track2()->Track()->Charge();
 193     double ptv1 = pair->Track1()->Track()->Pt();
 194     double ptv2 = pair->Track2()->Track()->Pt();
 195     double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
 196     double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
 197     double arg1 = -0.3 * 0.5 * chg1 * fRadius/(2*ptv1);
 198     double arg2 = -0.3 * 0.5 * chg2 * fRadius/(2*ptv2);
 199     double phid = phi2 - phi1 + TMath::ASin(arg2) - TMath::ASin(arg1);
 200
 201     while (phid>pih) phid -= pit;
 202     while (phid<-pih) phid += pit;
 203     //    dist = phi2 - phi1 + TMath::ASin(-0.3 * 0.5 * chg2 * fRadius/(2*ptv2)) - TMath::ASin(-0.3 * 0.5 * chg1 * fRadius/(2*ptv1));
 204     double etad = eta2 - eta1;
 205     fRadDDenominator->Fill(phid, etad);
 206   }
 207 }
 208
 209
 210 void AliFemtoTPCInnerCorrFctn::WriteHistos()
 211 {
 212   // Write out result histograms
 213   fDTPCNumerator->Write();
 214   fDTPCDenominator->Write();
 215   fRadDNumerator->Write();
 216   fRadDDenominator->Write();
 217 }
 218 //______________________________
 219 TList* AliFemtoTPCInnerCorrFctn::GetOutputList()
 220 {
 221   // Prepare the list of objects to be written to the output
 222   TList *tOutputList = new TList();
 223
 224   tOutputList->Add(fDTPCNumerator);
 225   tOutputList->Add(fDTPCDenominator);
 226   tOutputList->Add(fRadDNumerator);
 227   tOutputList->Add(fRadDDenominator);
 228
 229   return tOutputList;
 230 }
 231
 232 void AliFemtoTPCInnerCorrFctn::SetRadius(double rad)
 233 {
 234   fRadius = rad;
 235 }