1 ///////////////////////////////////////////////////////////////////////////
3 // AliFemtoQinvCorrFctn: //
4 // a simple Q-invariant correlation function //
6 ///////////////////////////////////////////////////////////////////////////
8 #include "AliFemtoQinvCorrFctn.h"
9 //#include "AliFemtoHisto.h"
13 ClassImp(AliFemtoQinvCorrFctn)
16 //____________________________
17 AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
22 fDetaDphiscal(kFALSE),
23 fPairKinematics(kFALSE),
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);
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);
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);
48 char tTitkT[101] = "kTDep";
49 strncat(tTitkT,title, 100);
50 fkTMonitor = new TH1D(tTitkT,title,250,0.0,5.0);
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);
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);
60 char tTitPair[101] = "Pair";
61 strncat(tTitPair,title, 100);
62 PairReader = new TNtuple(tTitPair,title, "px1:py1:pz1:e1:px2:py2:pz2:e2");
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);
71 // to enable error bar calculation...
73 fDenominator->Sumw2();
77 fNumDEtaDPhiS->Sumw2();
78 fDenDEtaDPhiS->Sumw2();
82 //____________________________
83 AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(const AliFemtoQinvCorrFctn& aCorrFctn) :
89 fDetaDphiscal(kFALSE),
90 fPairKinematics(kFALSE),
100 fNumerator = new TH1D(*aCorrFctn.fNumerator);
101 fDenominator = new TH1D(*aCorrFctn.fDenominator);
102 fRatio = new TH1D(*aCorrFctn.fRatio);
103 fkTMonitor = new TH1D(*aCorrFctn.fkTMonitor);
105 fNumDEtaDPhiS = new TH2D(*aCorrFctn.fNumDEtaDPhiS);
106 fDenDEtaDPhiS = new TH2D(*aCorrFctn.fDenDEtaDPhiS);
108 fDetaDphiscal = aCorrFctn.fDetaDphiscal;
109 fRaddedps = aCorrFctn.fRaddedps;
111 fPairKinematics = aCorrFctn.fPairKinematics;
113 if (aCorrFctn.PairReader)
114 PairReader = (TNtuple*)aCorrFctn.PairReader;
117 //____________________________
118 AliFemtoQinvCorrFctn::~AliFemtoQinvCorrFctn(){
124 delete fNumDEtaDPhiS;
125 delete fDenDEtaDPhiS;
129 //_________________________
130 AliFemtoQinvCorrFctn& AliFemtoQinvCorrFctn::operator=(const AliFemtoQinvCorrFctn& aCorrFctn)
132 // assignment operator
133 if (this == &aCorrFctn)
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);
145 if (fNumDEtaDPhiS) delete fNumDEtaDPhiS;
146 fNumDEtaDPhiS = new TH2D(*aCorrFctn.fNumDEtaDPhiS);
147 if (fDenDEtaDPhiS) delete fDenDEtaDPhiS;
148 fDenDEtaDPhiS = new TH2D(*aCorrFctn.fDenDEtaDPhiS);
150 fDetaDphiscal = aCorrFctn.fDetaDphiscal;
151 fRaddedps = aCorrFctn.fRaddedps;
153 fPairKinematics = aCorrFctn.fPairKinematics;
155 if (aCorrFctn.PairReader)
156 PairReader = (TNtuple*)aCorrFctn.PairReader;
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();
170 fRatio->Divide(fNumerator,fDenominator,1.0,1.0);
174 //____________________________
175 AliFemtoString AliFemtoQinvCorrFctn::Report(){
177 string stemp = "Qinv Correlation Function Report:\n";
179 snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
181 snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
183 snprintf(ctemp , 100, "Number of entries in ratio:\t%E\n",fRatio->GetEntries());
185 // stemp += mCoulombWeight->Report();
186 AliFemtoString returnThis = stemp;
189 //____________________________
190 void AliFemtoQinvCorrFctn::AddRealPair(AliFemtoPair* pair){
193 if (!fPairCut->Pass(pair)) return;
195 double tQinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs...
197 fNumerator->Fill(tQinv);
198 fkTMonitor->Fill(pair->KT());
201 //_______________________________________
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();
213 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
217 //AliWarning("Could not get AODInputHandler");
220 AliAODEvent *fAOD; // = new AliAODEvent()
221 fAOD = aodH->GetEvent();
222 magsign = fAOD->GetMagneticField();
230 else if ( magsign < 1)
235 Double_t rad = fRaddedps;
237 double afsi0b = 0.07510020733*chg1*fMagSign*rad/ptv1;
238 double afsi1b = 0.07510020733*chg2*fMagSign*rad/ptv2;
239 Double_t dps6 = phi2 - phi1 + TMath::ASin(afsi1b) - TMath::ASin(afsi0b);
240 dps6 = TVector2::Phi_mpi_pi(dps6);
242 // Double_t dps = (phi1-phi2+(TMath::ASin(-0.075*chg1*fMagSign*rad/ptv1))-(TMath::ASin(-0.075*chg2*fMagSign*rad/ptv2)));
243 // dps = TVector2::Phi_mpi_pi(dps);
244 double etad = eta2 - eta1;
246 fNumDEtaDPhiS->Fill(dps6,etad);
248 //_______________________________________________________________
250 // cout << "AliFemtoQinvCorrFctn::AddRealPair : " << pair->qInv() << " " << tQinv <<
251 //" " << pair->track1().FourMomentum() << " " << pair->track2().FourMomentum() << endl;
254 //____________________________
255 void AliFemtoQinvCorrFctn::AddMixedPair(AliFemtoPair* pair){
256 // add mixed (background) pair
258 if (!fPairCut->Pass(pair)) return;
261 double tQinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs...
262 fDenominator->Fill(tQinv,weight);
264 if (fPairKinematics) {
265 AliFemtoParticle* fTrack1 = pair->Track1();
266 AliFemtoParticle* fTrack2 = pair->Track2();
268 double px1 = fTrack1->FourMomentum().vect().x();
269 double py1 = fTrack1->FourMomentum().vect().y();
270 double pz1 = fTrack1->FourMomentum().vect().z();
271 double e1 = fTrack1->FourMomentum().e();
273 double px2 = fTrack2->FourMomentum().vect().x();
274 double py2 = fTrack2->FourMomentum().vect().y();
275 double pz2 = fTrack2->FourMomentum().vect().z();
276 double e2 = fTrack2->FourMomentum().e();
277 PairReader->Fill(px1, py1, pz1, e1, px2, py2, pz2, e2);
280 //_______________________________________
283 double phi1 = pair->Track1()->Track()->P().Phi();
284 double phi2 = pair->Track2()->Track()->P().Phi();
285 double chg1 = pair->Track1()->Track()->Charge();
286 double chg2 = pair->Track2()->Track()->Charge();
287 double ptv1 = pair->Track1()->Track()->Pt();
288 double ptv2 = pair->Track2()->Track()->Pt();
289 double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
290 double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
292 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
297 //AliWarning("Could not get AODInputHandler");
298 cout << "Could not get AODInputHandler" << endl;
302 fAOD = aodH->GetEvent();
303 magsign = fAOD->GetMagneticField();
310 else if ( magsign < 1)
315 Double_t rad = fRaddedps;
317 double afsi0b = 0.07510020733*chg1*fMagSign*rad/ptv1;
318 double afsi1b = 0.07510020733*chg2*fMagSign*rad/ptv2;
319 Double_t dps6 = phi2 - phi1 + TMath::ASin(afsi1b) - TMath::ASin(afsi0b);
320 dps6 = TVector2::Phi_mpi_pi(dps6);
321 double etad = eta2 - eta1;
323 fDenDEtaDPhiS->Fill(dps6,etad);
325 //_______________________________________________________________
328 //____________________________
329 void AliFemtoQinvCorrFctn::Write(){
330 // Write out neccessary objects
332 fDenominator->Write();
335 fNumDEtaDPhiS->Write();
336 fDenDEtaDPhiS->Write();
338 if (fPairKinematics) {
342 //______________________________
343 TList* AliFemtoQinvCorrFctn::GetOutputList()
345 // Prepare the list of objects to be written to the output
346 TList *tOutputList = new TList();
348 tOutputList->Add(fNumerator);
349 tOutputList->Add(fDenominator);
350 tOutputList->Add(fkTMonitor);
352 tOutputList->Add(fNumDEtaDPhiS);
353 tOutputList->Add(fDenDEtaDPhiS);
355 if (fPairKinematics) {
356 tOutputList->Add(PairReader);
361 void AliFemtoQinvCorrFctn::CalculateDetaDphis(Bool_t dedpsc, Double_t rad) {
362 fDetaDphiscal = dedpsc;
366 void AliFemtoQinvCorrFctn::CalculatePairKinematics(Bool_t pk) {
367 fPairKinematics = pk;