JCORRAN code update from DongJo
[u/mrichter/AliRoot.git] / PWGCF / Correlations / JCORRAN / AliJCorrelations.cxx
CommitLineData
66be7134 1/**************************************************************************
2 * Copyright(c) 1998-2014, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16// Comment describing what this class does needed!
17
18#include "AliJCorrelations.h"
19#include "AliJDataManager.h"
20
21#include "AliJCorrelations.h"
22#include "AliJHistos.h"
23#include "AliJBaseTrack.h"
24#include "AliJCard.h"
25#include "AliJHistos.h"
26#include "AliJRunTable.h"
27
28
29
30//==========================================================================
31// don't count the trigger here ! You'll miss the not associated triggers
32// but you will miss them in the main code too, unles you loop over those
33// fevent with no cgl (set the masMix to 1)
34//==========================================================================
35// blah
36
37AliJCorrelations::AliJCorrelations( AliJCard *cardIn, AliJHistos *histosIn) :
38 fcard(cardIn),
39 fhistos(histosIn),
40 fnReal(0),
41 fnMix(0),
42 fsumTriggerAndAssoc(0),
43 fsamplingMethod(0), // flat by default
44 fIsHeavyIon(0),
45 fawayPhiGap(0),
46 fmaxEtaRange(0),
47 fRSignalBin(0),
48 frandom(0x0),
49 fptt(0),
50 fpta(0),
51 fTrackPairEfficiency(0),
52 fIsIsolatedTrigger(false),
53 fpttBin(0),
54 fptaBin(0),
55 fPhiTrigger(0),
56 fPhiAssoc(0),
57 fDeltaPhi(0),
58 fDeltaPhiPiPi(0),
59 fDeltaEta(0),
60 fXlong(0),
61 fNearSide(true),
62 fEtaGapBin(0),
63 fPhiGapBinNear(0),
64 fPhiGapBinAway(0),
65 fRGapBinNear(0),
66 fRGapBinAway(0),
67 fCentralityBin(0),
68 fXlongBin(0),
69 fGeometricAcceptanceCorrection(1)
70{
71 // constructor
72
73 fsumTriggerAndAssoc = int( fcard->Get("sumTriggerAndAssoc") );
74 fmaxEtaRange = fcard->Get("EtaRange");
75 fRSignalBin = int(fcard->Get("EtaGapSignalBin"));
76
77 fDPhiUERegion[0] = fcard->Get("DPhiUERegion",0);
78 fDPhiUERegion[1] = fcard->Get("DPhiUERegion",1);
79 cout << fmaxEtaRange <<" fDPhiUERegion[0]="<< fDPhiUERegion[0] <<" fDPhiUERegion[1]="<< fDPhiUERegion[1] <<endl;
80 fIsHeavyIon = AliJRunTable::GetInstance().IsHeavyIon();
81
82 // -----------------------------------------------------------------------------------------------
83 // HARD CODED NUMBERS - VIOLATIONS - BREAKS the code when only on bin used in card.input!!!
84 // Who did this? Further more it us unused variable !!!!
85 // Jan
86 //fawayPhiGap = fcard->GetBinBorder(kEtaGapType, 2);
87 // -----------------------------------------------------------------------------------------------
88
89 for(int iRGap=0; iRGap < fcard->GetNoOfBins(kRGapType); iRGap++)
90 fRGap[iRGap] = fcard->GetBinBorder( kRGapType, iRGap);
91
92 //dPhiRange = fhistos->GetDphiRange();
93 frandom = new TRandom3(); //FK// frandom generator for jt flow UE
94 frandom->SetSeed(0); //FK//
95
96}
97
98AliJCorrelations::AliJCorrelations() :
99 fcard(0x0),
100 fhistos(0x0),
101 fnReal(0),
102 fnMix(0),
103 fsumTriggerAndAssoc(0),
104 fsamplingMethod(0), // flat by default
105 fIsHeavyIon(0),
106 fawayPhiGap(0),
107 fmaxEtaRange(0),
108 fRSignalBin(0),
109 frandom(0x0),
110 fptt(0),
111 fpta(0),
112 fTrackPairEfficiency(0),
113 fIsIsolatedTrigger(false),
114 fpttBin(0),
115 fptaBin(0),
116 fPhiTrigger(0),
117 fPhiAssoc(0),
118 fDeltaPhi(0),
119 fDeltaPhiPiPi(0),
120 fDeltaEta(0),
121 fXlong(0),
122 fNearSide(true),
123 fEtaGapBin(0),
124 fPhiGapBinNear(0),
125 fPhiGapBinAway(0),
126 fRGapBinNear(0),
127 fRGapBinAway(0),
128 fCentralityBin(0),
129 fXlongBin(0),
130 fGeometricAcceptanceCorrection(1)
131{
132 // default constructor
133}
134
135AliJCorrelations::AliJCorrelations(const AliJCorrelations& in) :
136 fcard(in.fcard),
137 fhistos(in.fhistos),
138 fnReal(in.fnReal),
139 fnMix(in.fnMix),
140 fsumTriggerAndAssoc(in.fsumTriggerAndAssoc),
141 fsamplingMethod(in.fsamplingMethod),
142 fIsHeavyIon(in.fIsHeavyIon),
143 fawayPhiGap(in.fawayPhiGap),
144 fmaxEtaRange(in.fmaxEtaRange),
145 fRSignalBin(in.fRSignalBin),
146 frandom(in.frandom),
147 fptt(in.fptt),
148 fpta(in.fpta),
149 fTrackPairEfficiency(in.fTrackPairEfficiency),
150 fIsIsolatedTrigger(in.fIsIsolatedTrigger),
151 fpttBin(in.fpttBin),
152 fptaBin(in.fptaBin),
153 fPhiTrigger(in.fPhiTrigger),
154 fPhiAssoc(in.fPhiAssoc),
155 fDeltaPhi(in.fDeltaPhi),
156 fDeltaPhiPiPi(in.fDeltaPhiPiPi),
157 fDeltaEta(in.fDeltaEta),
158 fXlong(in.fXlong),
159 fNearSide(in.fNearSide),
160 fEtaGapBin(in.fEtaGapBin),
161 fPhiGapBinNear(in.fPhiGapBinNear),
162 fPhiGapBinAway(in.fPhiGapBinAway),
163 fRGapBinNear(in.fRGapBinNear),
164 fRGapBinAway(in.fRGapBinAway),
165 fCentralityBin(in.fCentralityBin),
166 fXlongBin(in.fXlongBin),
167 fGeometricAcceptanceCorrection(in.fGeometricAcceptanceCorrection)
168{
169 // The pointers to card and histos are just copied. I think this is safe, since they are not created by
170 // AliJCorrelations and thus should not disappear if the AliJCorrelation managing them is destroyed.
171
172 fDPhiUERegion[0] = in.fDPhiUERegion[0];
173 fDPhiUERegion[1] = in.fDPhiUERegion[1];
174
175 for(int iRGap=0; iRGap < fcard->GetNoOfBins(kRGapType); iRGap++){
176 fRGap[iRGap] = in.fRGap[iRGap];
177 }
178
179 frandom = new TRandom3(); // frandom generator for jt flow UE
180 frandom->SetSeed(0);
181}
182
183AliJCorrelations& AliJCorrelations::operator=(const AliJCorrelations& in){
184 // Assingment operator
185
186 if (&in==this) return *this;
187
188 fptt = in.fptt;
189 fpta = in.fpta;
190 fTrackPairEfficiency = in.fTrackPairEfficiency;
191 fIsIsolatedTrigger = in.fIsIsolatedTrigger;
192 fpttBin = in.fpttBin;
193 fptaBin = in.fptaBin;
194 fPhiTrigger = in.fPhiTrigger;
195 fPhiAssoc = in.fPhiAssoc;
196 fDeltaPhi = in.fDeltaPhi;
197 fDeltaPhiPiPi = in.fDeltaPhiPiPi;
198 fDeltaEta = in.fDeltaEta;
199 fXlong = in.fXlong;
200 fNearSide = in.fNearSide;
201 fEtaGapBin = in.fEtaGapBin;
202 fPhiGapBinNear = in.fPhiGapBinNear;
203 fPhiGapBinAway = in.fPhiGapBinAway;
204 fRGapBinNear = in.fRGapBinNear;
205 fRGapBinAway = in.fRGapBinAway;
206 fCentralityBin = in.fCentralityBin;
207 fXlongBin = in.fXlongBin;
208 fGeometricAcceptanceCorrection = in.fGeometricAcceptanceCorrection;
209 fnReal = in.fnReal;
210 fnMix = in.fnMix;
211 fsumTriggerAndAssoc = in.fsumTriggerAndAssoc;
212 fsamplingMethod = in.fsamplingMethod;
213 fIsHeavyIon = in.fIsHeavyIon;
214 fawayPhiGap = in.fawayPhiGap;
215 fmaxEtaRange = in.fmaxEtaRange;
216 fRSignalBin = in.fRSignalBin;
217
218 // The pointers to card and histos are just copied. I think this is safe, since they are not created by
219 // AliJCorrelations and thus should not disappear if the AliJCorrelation managing them is destroyed.
220 fcard = in.fcard;
221 fhistos = in.fhistos;
222
223 fDPhiUERegion[0] = in.fDPhiUERegion[0];
224 fDPhiUERegion[1] = in.fDPhiUERegion[1];
225
226 for(int iRGap=0; iRGap < fcard->GetNoOfBins(kRGapType); iRGap++){
227 fRGap[iRGap] = in.fRGap[iRGap];
228 }
229
230 frandom = new TRandom3(); // frandom generator for jt flow UE
231 frandom->SetSeed(0);
232
233 return *this;
234 // copy constructor
235}
236
237
238void AliJCorrelations::FillHisto(corrFillType cFTyp, fillType fTyp, int cBin, int zBin, AliJBaseTrack *ftk1, AliJBaseTrack *ftk2){
239 // histo filler
240 if( cFTyp == kAzimuthFill )
241 FillAzimuthHistos( fTyp, cBin, zBin, ftk1, ftk2);
242
243}
244
245//=============================================================================================
246void AliJCorrelations::FillAzimuthHistos(fillType fTyp, int CentBin, int ZBin, AliJBaseTrack *ftk1, AliJBaseTrack *ftk2)
247//=============================================================================================
248{
249 // histo filler
250 bool twoTracks = false;
251 if(ftk1->GetParticleType()==kHadron && ftk2->GetParticleType()==kHadron) twoTracks =true;
252
253 //double-counting check
254 if(fTyp == kReal && twoTracks && ftk1->GetID()==ftk2->GetID()) return;
255
256 // Accept only unlike sign pairs
257 //if(ftk1->GetCharge() > 0 && ftk2->GetCharge() > 0) return;
258 //if(ftk1->GetCharge() < 0 && ftk2->GetCharge() < 0) return;
259
260 //----------------------------------------------------------------
261 fptt = ftk1->Pt();
262 fpta = ftk2->Pt();
263
264 fTrackPairEfficiency = 1./( ftk1->GetTrackEff() * ftk2->GetTrackEff() );
265
266 fIsIsolatedTrigger = ftk1->GetIsIsolated()>0 ? true : false; //FK// trigger particle is isolated hadron
267
268 //phit= ftk1->Phi(); //for RP
269 fpttBin = ftk1->GetTriggBin();
270 fptaBin = ftk2->GetAssocBin();
271 fPhiTrigger = ftk1->Phi();
272 fPhiAssoc = ftk2->Phi();
273 fDeltaPhi = DeltaPhi(fPhiTrigger, fPhiAssoc); //radians
274 fDeltaPhiPiPi = atan2(sin(fPhiTrigger-fPhiAssoc), cos(fPhiTrigger-fPhiAssoc));
275 fDeltaEta = ftk1->Eta() - ftk2->Eta();
276 //double dEtaFar = ftk1->Eta() + ftk2->Eta();
277
278 fNearSide = cos(fPhiTrigger-fPhiAssoc) > 0 ? true : false;
279
280 fEtaGapBin = fcard->GetBin( kEtaGapType, fabs(fDeltaEta));
281 fPhiGapBinNear = fcard->GetBin( kEtaGapType, fabs(fDeltaPhiPiPi) );
282 fPhiGapBinAway = fcard->GetBin( kEtaGapType, fabs(fDeltaPhi-kJPi) ); //here the angle must be 0-2pi and not (-pi,pi)
283 fRGapBinNear = fcard->GetBin( kRGapType, fabs(ftk1->DeltaR(*ftk2)));
284 fRGapBinAway = fcard->GetBin( kRGapType, sqrt(pow(fDeltaPhi-kJPi,2)+fDeltaEta*fDeltaEta) );
285 fCentralityBin = CentBin;
286
287 TLorentzVector vTrigger = ftk1->GetLorentzVector(), vAssoc = ftk2->GetLorentzVector(); // Lorentz vectors for trigger and associated particles
288 fXlong = vTrigger.Vect().Dot(vAssoc.Vect())/pow(vTrigger.P(),2);
289 fXlongBin = fcard->GetBin(kXeType, TMath::Abs(fXlong));
290
291 //if( rGapBin != fRGapBinNear ) cout<<"dR vs fRGapBinNear = "<<rGapBin<<"\t"<<fRGapBinNear<<endl;
292 //if( rGapBin != fRGapBinAway ) cout<<"dR vs fRGapBinAway = "<<rGapBin<<"\t"<<fRGapBinAway<<endl;
293 //----------------------------------------------------------------
294
295 //acceptance correction triangle or mixed fevent
296 // fGeometricAcceptanceCorrection = 1;
297 fGeometricAcceptanceCorrection = ( fsamplingMethod == 0 ) ? GetGeoAccCorrFlat(fDeltaEta) : GetGeoAccCorrIncl(fDeltaEta);
298
299
300 if(fpttBin<0 || fptaBin<0 || fEtaGapBin<0 ){
301 cout<<"Error in FillAzimuthHistos: some pT or eta out of bin. pttBin="<<fpttBin<<" pTaBin="<<fptaBin <<" etaGapBin="<< fEtaGapBin << endl;
302 ftk1->Print();
303 ftk2->Print();
304 exit(-1);
305 //return;
306 }
307
308 if(fDeltaPhi==0) cout <<" fdphi=0; fptt="<< fptt<<" fpta="<<fpta<<" TID="<<ftk1->GetID()<<" AID="<<ftk2->GetID() <<" tphi="<< fPhiTrigger <<" aphi="<< fPhiAssoc << endl;
309
310 // ===================================================================
311 // ===================== Fill Histograms ===========================
312 // ===================================================================
313
314 bool fill2DBackgroundQualityControlHistograms = false; // Choose whether to fill the DeltaPhi DeltaEta histograms for jT background
315
316 //FillPairPtAndCosThetaStarHistograms(fTyp, ftk1, ftk2); // Fill the pair pT and cos(theta*) histograms TODO: Does not work! Needs debugging
317 FillXeHistograms(fTyp); // Fill the xE and xLong histograms
318 FillJtHistograms(fTyp, ftk1, ftk2, fill2DBackgroundQualityControlHistograms); // Fill the jT and pout histograms togerher with some background quality assurance histograms
319 FillDeltaEtaHistograms(fTyp, ZBin); // Fill all the delta eta histograms
320 FillDeltaPhiHistograms(fTyp); // Fill the azimuthal correlation functions
321 FillPtaHistograms(fTyp); // Fill various pTa histograms
322 FillIAAAndMoonHistograms(fTyp, ZBin); // Fill the I_AA and moon histograms
323
324}
325
326
327void AliJCorrelations::FillPairPtAndCosThetaStarHistograms(fillType fTyp, AliJBaseTrack *ftk1, AliJBaseTrack *ftk2)
328{
329 // This method fills the pair pT and Cos(thata*) histograms
330
331 TVector3 v3trigg = ftk1->GetLorentzVector().Vect();
332 TVector3 v3assoc = ftk2->GetLorentzVector().Vect();
333 double pairMass = sqrt(2*v3trigg.Mag()*v3assoc.Mag()*(1-cos(v3trigg.Angle(v3assoc))));
334 int imass = fcard->GetBin(kMassType, pairMass);
335 TVector3 v3pairPt = v3trigg + v3assoc;
336 double pairPt = v3pairPt.Perp();
337 double wPairPt = pairPt>0 ? 1/pairPt : 0;
338
339 //================================
340 // Pair pT and cos(Theta^star) === // TODO: Needs debugging, the code below does not work!
341 //================================
342
343 if(!fNearSide) {
344 fhistos->fhPairPt[fTyp][fpttBin][fptaBin]->Fill( pairPt, fTrackPairEfficiency);
345 double pairDphi = v3pairPt.DeltaPhi(v3trigg);
346 if(pairDphi<-kJPi/3.0) pairDphi += kJTwoPi;
347 fhistos->fhPairPtDphi[fTyp][fpttBin][fptaBin]->Fill( pairDphi/kJPi );
348 }
349
350 if ( fTyp == kReal ) {
351 //cout<<" fptt="<< fptt <<" fpta="<< fpta <<" fdphi="<< fDeltaPhi <<" pairPt="<< pairPt << endl;
352
353 int ipairPt = fcard->IsLessThanUpperPairPtCut(pairPt);
354 //cout<<"ppt="<<pairPt<<" mass="<<pairMass<<" ip="<<ipairPt<<" im="<<imass<<endl;
355 if(imass>=0) {
356 fhistos->fhPairPtMass[imass]->Fill( pairPt, fTrackPairEfficiency*wPairPt );
357 fhistos->fhPairDPhi[imass]->Fill( fDeltaPhi/kJPi, fTrackPairEfficiency );
358 fhistos->fhPairDpT[imass]->Fill( fptt-fpta, fTrackPairEfficiency );
359 }
360 if(imass>=0 && ipairPt>=0 ){
361 double cmsRap = (ftk1->Eta()+ftk2->Eta())/2.;
362 double cosThetaStar = fabs(cos(2*atan(exp(-ftk1->Eta()+cmsRap))));
363 fhistos->fhCosThetaStar[fTyp][ipairPt][imass]->Fill( cosThetaStar, fTrackPairEfficiency );
364
365 //frandom boost, Let's try just trigger phi as a cms boost.
366 cosThetaStar = fabs(cos(2*atan(exp(-ftk1->Eta()+fPhiTrigger))));
367 fhistos->fhCosThetaStar[kMixed][ipairPt][imass]->Fill( cosThetaStar, fTrackPairEfficiency );
368
369 fhistos->fhInvMass[ipairPt]->Fill(pairMass, fTrackPairEfficiency);
370 fhistos->fhCMSrap[ipairPt][imass]->Fill( cosThetaStar, cmsRap*fTrackPairEfficiency );
371 fhistos->fpCMSrap->Fill( cosThetaStar, fabs(cmsRap)*fTrackPairEfficiency );
372 }
373 }
374
375 //---------------------------------
376 //if(pairMass>10 && fcard->IsLessThanUpperPairPtCut(pairPt) ){
377 // cout<<pairMass<<" e1="<<ftk1->Eta()-cmsRap<<" e2="<<ftk2->Eta()-cmsRap<<" cms="<<cmsRap;
378 // cout<<" ppt="<<pairPt<<" tphi="<<fPhiTrigger<<" aphi="<<fPhiAssoc<<" t1="<<fptt<<" t2="<<fpta<<endl;
379 // double thetaStar = fabs(2*atan(exp(-ftk1->Eta()-cmsRap)));
380 // cout<<thetaStar<<" "<<fabs(cos(thetaStar))<<endl;
381 //}
382
383 //---------------------------------
384}
385
386void AliJCorrelations::FillXeHistograms(fillType fTyp)
387{
388 // This method fills the xE and xLong histograms
389
390 double xe = -fpta*cos(fPhiTrigger-fPhiAssoc)/fptt;
391
392 if( fTyp == kReal ) {
393 fhistos->fhxEPtBin[0][fpttBin][fptaBin]->Fill(fXlong, fGeometricAcceptanceCorrection * fTrackPairEfficiency);
394 if( fNearSide ) {
395 fhistos->fhxEPtBin[1][fpttBin][fptaBin]->Fill(fXlong, fGeometricAcceptanceCorrection * fTrackPairEfficiency);
396 } else {
397 fhistos->fhxEPtBin[2][fpttBin][fptaBin]->Fill(fXlong, fGeometricAcceptanceCorrection * fTrackPairEfficiency);
398 }
399 }
400
401 if(fNearSide) {
402 fhistos->fhxEN [fTyp][fpttBin]->Fill(-xe, fGeometricAcceptanceCorrection * fTrackPairEfficiency);
403 } else {
404 fhistos->fhxEF [fTyp][fpttBin]->Fill(xe, fGeometricAcceptanceCorrection * fTrackPairEfficiency);
405 if(fIsIsolatedTrigger) fhistos->fhxEFIsolTrigg[fTyp][fpttBin]->Fill(xe, fGeometricAcceptanceCorrection * fTrackPairEfficiency);
406 }
407}
408
409void AliJCorrelations::FillJtHistograms(fillType fTyp, AliJBaseTrack *ftk1, AliJBaseTrack *ftk2, bool fill2DBackground)
410{
411 // This method fills the jT and pout histograms
412
413 double pout = fabs(fpta*sin(fPhiTrigger-fPhiAssoc));
414
415 if(fNearSide) {
416
417 TLorentzVector vTrigger = ftk1->GetLorentzVector(), vAssoc = ftk2->GetLorentzVector(); // Lorentz vectors for tracks
418
419 double jt = vAssoc.Perp(vTrigger.Vect());
420 double weight = jt > 1e-3 ? fGeometricAcceptanceCorrection * fTrackPairEfficiency/jt : 0;
421
422 int iKlongBin = fcard->GetBin(kLongType, vTrigger.Vect().Dot(vAssoc.Vect())/vTrigger.P());
423 if(iKlongBin>=0){
424 fhistos->fhJTKlong[fTyp][fCentralityBin][fpttBin][iKlongBin]->Fill(jt, weight);
425 //cout <<" jt="<< jt <<" geo="<< fGeometricAcceptanceCorrection <<" eff="<< fTrackPairEfficiency <<" w="<< fGeometricAcceptanceCorrection * fTrackPairEfficiency/jt <<endl;
426 }
427 if(fXlongBin>=0){
428 fhistos->fhJT[fTyp][fCentralityBin][fpttBin][fXlongBin]->Fill(jt, weight);
429 //cout <<" jt="<< jt <<" geo="<< fGeometricAcceptanceCorrection <<" eff="<< fTrackPairEfficiency <<" w="<< fGeometricAcceptanceCorrection * fTrackPairEfficiency/jt <<endl;
430 }
431 if(fptaBin>=0){
432 fhistos->fhJTPta[fTyp][fCentralityBin][fpttBin][fptaBin]->Fill(jt, weight);
433 //cout <<" jt="<< jt <<" geo="<< fGeometricAcceptanceCorrection <<" eff="<< fTrackPairEfficiency <<" w="<< fGeometricAcceptanceCorrection * fTrackPairEfficiency/jt <<endl;
434 }
435
436 TLorentzVector vAssocRndm;
437 TLorentzVector vThrustRndm;
438
439 //-------------------------------------
440 // Rapidity Gap background
441 // ------------------------------------
442 double etaTriggRndm, etaAssocRndm;
443 double dphiAssocRndm = 0;
444
445 if(fTyp==kReal && ( fEtaGapBin >=0 || fRGapBinNear >=0 ) ){
446 for(int iEtaGap=0; iEtaGap<=fEtaGapBin; iEtaGap++) fhistos->fhPtaEtaGapN[fCentralityBin][iEtaGap][fpttBin]->Fill(fpta, fTrackPairEfficiency);
447 for(int iRGap=0; iRGap<=fRGapBinNear; iRGap++) fhistos->fhPtaRGapN[fCentralityBin][iRGap][fpttBin]->Fill(fpta, fTrackPairEfficiency);
448 for(int itrial=0; itrial<20; itrial++){ //For each high eta gap track generate ten others
449
450 if(fsamplingMethod == 0){
451 etaTriggRndm = frandom->Uniform(-fmaxEtaRange, fmaxEtaRange);
452 etaAssocRndm = frandom->Uniform(-fmaxEtaRange, fmaxEtaRange);
453 } else {
454 etaTriggRndm = fhistos->fhIetaTriggFromFile[fCentralityBin][fpttBin]->GetRandom();
455 etaAssocRndm = fhistos->fhIetaAssocFromFile[fCentralityBin][fptaBin]->GetRandom();
456 }
457 if( fEtaGapBin >=0 || fRGapBinNear >=0 ) {
458 if( fsamplingMethod == 0) {
459 dphiAssocRndm = kJPi * frandom->Uniform(-0.5, 0.5);
460 } else {
461 dphiAssocRndm = fhistos->fhIphiAssocFromFile[fCentralityBin][fptaBin]->GetRandom();
462 }
463 vAssocRndm.SetPtEtaPhiM(fpta, etaAssocRndm, fPhiTrigger + dphiAssocRndm, 0); //set randomized assoc TLorentz vector
464 } else {
465 vAssocRndm.SetPtEtaPhiM(fpta, etaAssocRndm, fPhiAssoc, 0); //set randomized assoc TLorentz vector
466 }
467
468 vThrustRndm.SetPtEtaPhiM(vTrigger.Pt(),etaTriggRndm, fPhiTrigger, 0);
469
470 double dEtaRndm = etaTriggRndm - etaAssocRndm;
471 double geoAccCorrRndm = (fsamplingMethod == 0 || fPhiGapBinNear<0 ) ? GetGeoAccCorrFlat(dEtaRndm) : GetGeoAccCorrIncl(dEtaRndm);
472
473 jt = vAssocRndm.Perp(vThrustRndm.Vect());
474
475 // TODO : shadowing of iKlongBin
476 int jKlongBin = fcard->GetBin(kLongType, vThrustRndm.Vect().Dot(vAssocRndm.Vect())/vThrustRndm.P());
477 if(jKlongBin>=0){ //
478 if(jt>1e-3) { //here we used and BgFidCut
479 for(int iEtaGap=0; iEtaGap<=fEtaGapBin; iEtaGap++){
480 fhistos->fhJTKlongBg[fCentralityBin][iEtaGap][fpttBin][jKlongBin]->Fill(jt, geoAccCorrRndm * fTrackPairEfficiency/jt);
481 if(fill2DBackground) fhistos->fhDphiDetaKlong[fCentralityBin][iEtaGap][fpttBin][jKlongBin]->Fill(dEtaRndm, dphiAssocRndm, geoAccCorrRndm * fTrackPairEfficiency);
482 fhistos->fhBgAssocKlong[fCentralityBin][iEtaGap][fpttBin][jKlongBin]->Fill(fpta, fTrackPairEfficiency);
483 }
484 for(int iRGap=0; iRGap<=fRGapBinNear; iRGap++){
485 fhistos->fhJTKlongBgR[fCentralityBin][iRGap][fpttBin][jKlongBin]->Fill(jt, geoAccCorrRndm * fTrackPairEfficiency/jt);
486 if(fill2DBackground) fhistos->fhDphiDetaKlongR[fCentralityBin][iRGap][fpttBin][jKlongBin]->Fill(dEtaRndm, dphiAssocRndm, geoAccCorrRndm * fTrackPairEfficiency);
487 fhistos->fhBgAssocKlongR[fCentralityBin][iRGap][fpttBin][jKlongBin]->Fill(fpta, fTrackPairEfficiency);
488 }
489 }
490 }
491 int iXeBin = fcard->GetBin(kXeType, vThrustRndm.Vect().Dot(vAssocRndm.Vect())/pow(vThrustRndm.P(),2) );
492 if(iXeBin>=0){
493 if(jt>1e-3) { //here we used and BgFidCut
494 for(int iEtaGap=0; iEtaGap<=fEtaGapBin; iEtaGap++){
495 fhistos->fhJTBg[fCentralityBin][iEtaGap][fpttBin][iXeBin]->Fill(jt, geoAccCorrRndm * fTrackPairEfficiency/jt);
496 if(fill2DBackground) fhistos->fhDphiDetaXe[fCentralityBin][iEtaGap][fpttBin][iXeBin]->Fill(dEtaRndm, dphiAssocRndm, geoAccCorrRndm * fTrackPairEfficiency);
497 fhistos->fhBgAssocXe[fCentralityBin][iEtaGap][fpttBin][iXeBin]->Fill(fpta, fTrackPairEfficiency);
498 }
499 for(int iRGap=0; iRGap<=fRGapBinNear; iRGap++){
500 fhistos->fhJTBgR[fCentralityBin][iRGap][fpttBin][iXeBin]->Fill(jt, geoAccCorrRndm * fTrackPairEfficiency/jt);
501 if(fill2DBackground) fhistos->fhDphiDetaXeR[fCentralityBin][iRGap][fpttBin][iXeBin]->Fill(dEtaRndm, dphiAssocRndm, geoAccCorrRndm * fTrackPairEfficiency);
502 fhistos->fhBgAssocXeR[fCentralityBin][iRGap][fpttBin][iXeBin]->Fill(fpta, fTrackPairEfficiency);
503 }
504 }
505 }
506 if(fptaBin>=0){
507 if(jt>1e-3) { //here we used and BgFidCut
508 for(int iEtaGap=0; iEtaGap<=fEtaGapBin; iEtaGap++){
509 fhistos->fhJTPtaBg[fCentralityBin][iEtaGap][fpttBin][fptaBin]->Fill(jt, geoAccCorrRndm * fTrackPairEfficiency/jt);
510 if(fill2DBackground) fhistos->fhDphiDetaPta[fCentralityBin][iEtaGap][fpttBin][fptaBin]->Fill(dEtaRndm, dphiAssocRndm, geoAccCorrRndm * fTrackPairEfficiency);
511 fhistos->fhBgAssocPta[fCentralityBin][iEtaGap][fpttBin][fptaBin]->Fill(fpta, fTrackPairEfficiency);
512 }
513 for(int iRGap=0; iRGap<=fRGapBinNear; iRGap++){
514 fhistos->fhJTPtaBgR[fCentralityBin][iRGap][fpttBin][fptaBin]->Fill(jt, geoAccCorrRndm * fTrackPairEfficiency/jt);
515 if(fill2DBackground) fhistos->fhDphiDetaPtaR[fCentralityBin][iRGap][fpttBin][fptaBin]->Fill(dEtaRndm, dphiAssocRndm, geoAccCorrRndm * fTrackPairEfficiency);
516 fhistos->fhBgAssocPtaR[fCentralityBin][iRGap][fpttBin][fptaBin]->Fill(fpta, fTrackPairEfficiency);
517 }
518 }
519 }
520 } //trials
521 }// Eta Gap Random
522
523
524 } else {
525 fhistos->fhPoutF[fTyp][fCentralityBin][fpttBin][fptaBin]->Fill(pout, fGeometricAcceptanceCorrection * fTrackPairEfficiency);
526 }
527}
528
529void AliJCorrelations::FillDeltaEtaHistograms(fillType fTyp, int ZBin)
530{
531 // This method fills the DeltaEta histograms
532
533 if( fNearSide ){ //one could check the phiGapBin, but in the pi/2 <1.6 and thus phiGap is always>-1
534 if( fTyp == 0 ) {
535 fhistos->fhDEtaNear[fCentralityBin][ZBin][fPhiGapBinNear][fpttBin][fptaBin]->Fill( fDeltaEta , fGeometricAcceptanceCorrection * fTrackPairEfficiency );
536 if(fXlongBin>=0) fhistos->fhDEtaNearXEbin[fCentralityBin][ZBin][fPhiGapBinNear][fpttBin][fXlongBin]->Fill( fDeltaEta , fGeometricAcceptanceCorrection * fTrackPairEfficiency );
537 } else {
538 fhistos->fhDEtaNearM[fCentralityBin][ZBin][fPhiGapBinNear][fpttBin][fptaBin]->Fill( fDeltaEta , fGeometricAcceptanceCorrection * fTrackPairEfficiency );
539 if(fXlongBin>=0) fhistos->fhDEtaNearMXEbin[fCentralityBin][ZBin][fPhiGapBinNear][fpttBin][fXlongBin]->Fill( fDeltaEta , fGeometricAcceptanceCorrection * fTrackPairEfficiency );
540 }
541 } else {
542 if(fPhiGapBinAway<=3) fhistos->fhDEtaFar[fTyp][fCentralityBin][fpttBin]->Fill( fDeltaEta, fGeometricAcceptanceCorrection * fTrackPairEfficiency );
543 }
544}
545
546void AliJCorrelations::FillDeltaPhiHistograms(fillType fTyp)
547{
548 // This method fills the DeltaPhi histograms
549
550 // When hists are filled for thresholds they are not properly normalized and need to be subtracted
551 // This induced improper errors - subtraction of not-independent entries
552
553 fhistos->fhDphiAssoc[fTyp][fCentralityBin][fEtaGapBin][fpttBin][fptaBin]->Fill( fDeltaPhi/kJPi , fGeometricAcceptanceCorrection * fTrackPairEfficiency);
554 if(fXlongBin>=0) fhistos->fhDphiAssocXEbin[fTyp][fCentralityBin][fEtaGapBin][fpttBin][fXlongBin]->Fill( fDeltaPhi/kJPi , fGeometricAcceptanceCorrection * fTrackPairEfficiency);
555
556 if(fIsIsolatedTrigger) fhistos->fhDphiAssocIsolTrigg[fTyp][fCentralityBin][fpttBin][fptaBin]->Fill( fDeltaPhi/kJPi , fGeometricAcceptanceCorrection * fTrackPairEfficiency); //FK//
557}
558
559void AliJCorrelations::FillPtaHistograms(fillType fTyp)
560{
561 // This method fills various pta histograms
562
563 if ( fTyp == kReal ) {
564 //must be here, not in main, to avoid counting triggers
565 fhistos->fhAssocPtBin[fCentralityBin][fpttBin][fptaBin]->Fill(fpta ); //I think It should not be weighted by Eff
566
567 //++++++++++++++++++++++++++++++++++++++++++++++++++
568 // in order to get mean pTa in the jet peak one has
569 // to fill fhMeanPtAssoc in |DeltaEta|<0.4
570 // +++++++++++++++++++++++++++++++++++++++++++++++++
571 if(fEtaGapBin>=0 && fEtaGapBin<2){
572 fhistos->fhMeanPtAssoc[fCentralityBin][fpttBin][fptaBin]->Fill( fDeltaPhi/kJPi , fpta );
573 fhistos->fhMeanZtAssoc[fCentralityBin][fpttBin][fptaBin]->Fill( fDeltaPhi/kJPi , fpta/fptt);
574 }
575
576 //UE distribution
577 if(fabs(fDeltaPhiPiPi/kJPi)>fDPhiUERegion[0] && fabs(fDeltaPhiPiPi/kJPi)<fDPhiUERegion[1]){
578 for(int iEtaGap=0; iEtaGap<=fEtaGapBin; iEtaGap++) //FK// UE Pta spectrum for different eta gaps
579 fhistos->fhPtAssocUE[fCentralityBin][iEtaGap][fpttBin]->Fill(fpta, fTrackPairEfficiency);
580 if(fIsIsolatedTrigger){ //FK// trigger is isolated hadron
581 fhistos->fhPtAssocUEIsolTrigg[fpttBin]->Fill(fpta, fTrackPairEfficiency); //FK//
582 }
583 }
584 if(fabs(fDeltaPhi/kJPi)<0.15) fhistos->fhPtAssocN[fpttBin]->Fill(fpta, fTrackPairEfficiency);
585 if(fabs(fDeltaPhi/kJPi-1)<0.15) fhistos->fhPtAssocF[fpttBin]->Fill(fpta, fTrackPairEfficiency);
586
587 fnReal++;
588 } else { // only mix
589 fnMix++;
590 }
591}
592
593void AliJCorrelations::FillIAAAndMoonHistograms(fillType fTyp, int ZBin)
594{
595 // This method fills the I_AA and moon histograms
596
597 fhistos->fhDphiAssoc2DIAA[fTyp][fCentralityBin][ZBin][fpttBin][fptaBin]->Fill( fDeltaEta, fDeltaPhi/kJPi, fTrackPairEfficiency);
598
599 if(fRGapBinNear>=0){
600 if(fRGapBinNear <= fRSignalBin) fhistos->fhDRNearPt[fTyp][fCentralityBin][ZBin][fRGapBinNear][fpttBin]->Fill( fpta, fGeometricAcceptanceCorrection * fTrackPairEfficiency );
601 // - moon -
602 if(fRGapBinNear>0){
603 for( int irs = 0; irs <= fRSignalBin;irs++ ){
604 if( fRGapBinNear < irs ) continue;
605 if( fPhiGapBinNear > irs ) continue;
606 for ( int ir1=0;ir1<= fRGapBinNear;ir1++ ){
607 if( irs > ir1 ) continue;
608 if( ir1 > fRGapBinNear ) continue;
609 double rGapS= fRGap[irs+1];
610 double rGap1= fRGap[ir1+1];
611 if( rGap1 < TMath::Abs(fDeltaPhiPiPi) ) continue;
612 if( rGapS < TMath::Abs(fDeltaPhiPiPi) ) continue;
613 double dEtaMin = sqrt(rGap1*rGap1-fDeltaPhiPiPi*fDeltaPhiPiPi);
614 double dEtaMax = dEtaMin + sqrt(rGapS*rGapS-fDeltaPhiPiPi*fDeltaPhiPiPi);
615 double eta = TMath::Abs( fDeltaEta );
616 if( eta > dEtaMin && eta < dEtaMax ){
617 //if( eta > dEtaMin && eta < dEtaMax && fDeltaEta <0 ){
618 // xxx
619 // fhistos->hDRNearPtMoon[fTyp][fCentralityBin][ZBin][ir1][irs][fpttBin]->Fill( fpta, fGeometricAcceptanceCorrection * fTrackPairEfficiency );
620
621 if( fTyp == 0 )
622 fhistos->fhDRNearPtMoon[fCentralityBin][ZBin][ir1][irs][fpttBin]->Fill( fpta, fGeometricAcceptanceCorrection * fTrackPairEfficiency );
623 else
624 fhistos->fhDRNearPtMoonM[fCentralityBin][ZBin][ir1][irs][fpttBin]->Fill( fpta, fGeometricAcceptanceCorrection * fTrackPairEfficiency );
625 if(fTyp == kReal) fhistos->fhDphiAssoc2D[ir1][irs]->Fill( fDeltaEta, fDeltaPhi/kJPi, fGeometricAcceptanceCorrection * fTrackPairEfficiency );
626 }
627 }
628 }
629 }
630 }
631
632 if(fRGapBinAway>=0){
633 if(fRGapBinAway <= fRSignalBin) fhistos->fhDRFarPt[fTyp][fCentralityBin][ZBin][fRGapBinAway][fpttBin]->Fill( fpta, fGeometricAcceptanceCorrection * fTrackPairEfficiency );
634 // - moon -
635 if(fRGapBinAway>0){
636 for( int irs = 0; irs <= fRSignalBin;irs++ ){
637 if( fRGapBinAway < irs ) continue;
638 if( fPhiGapBinAway > irs ) continue;
639 for ( int ir1=0;ir1<= fRGapBinAway;ir1++ ){
640 if( irs > ir1 ) continue;
641 if( ir1 > fRGapBinAway ) continue;
642 double rGapS= fRGap[irs+1];
643 double rGap1= fRGap[ir1+1];
644 if( rGap1 < TMath::Abs(fDeltaPhi-kJPi) ) continue;
645 if( rGapS < TMath::Abs(fDeltaPhi-kJPi) ) continue;
646 double dEtaMin = sqrt(rGap1*rGap1- (fDeltaPhi-kJPi)*(fDeltaPhi-kJPi) );
647 double dEtaMax = dEtaMin + sqrt(rGapS*rGapS-(fDeltaPhi-kJPi)*(fDeltaPhi-kJPi) );
648 double eta = TMath::Abs( fDeltaEta );
649 if( eta > dEtaMin && eta < dEtaMax ){
650 //if( eta > dEtaMin && eta < dEtaMax && fDeltaEta <0 ){
651 // xxx
652 // fhistos->hDRFarPtMoon[fTyp][fCentralityBin][ZBin][ir1][irs][fpttBin]->Fill( fpta, fGeometricAcceptanceCorrection * fTrackPairEfficiency );
653 if( fTyp == 0 )
654 fhistos->fhDRFarPtMoon[fCentralityBin][ZBin][ir1][irs][fpttBin]->Fill( fpta, fGeometricAcceptanceCorrection * fTrackPairEfficiency );
655 else
656 fhistos->fhDRFarPtMoonM[fCentralityBin][ZBin][ir1][irs][fpttBin]->Fill( fpta, fGeometricAcceptanceCorrection * fTrackPairEfficiency );
657
658 if(fTyp == kReal) fhistos->fhDphiAssoc2D[ir1][irs]->Fill( fDeltaEta, fDeltaPhi/kJPi, fGeometricAcceptanceCorrection * fTrackPairEfficiency );
659 }
660 }
661 }
662 }
663 }
664}
665
666
667double AliJCorrelations::GetGeoAccCorrFlat(double deltaEta){
668 //FK// calculate acceptance correction on pseudorapidity triangle
669 double absDEta = fabs(deltaEta);
670
671 double denominator = 1 - absDEta/(2*fmaxEtaRange);
672 //double denominator = 1 - (absDEta - ftriggFiducCut)/(2*fmaxEtaRange-2*ftriggFiducCut);//When Fid>0 max_Deta je mensi nez 2*EtaMax
673 //double denominator = 1 - (absDEta - ftriggFiducCut)/(2*fmaxEtaRange-ftriggFiducCut);
674 if(denominator > 1e-6)
675 return 1.0/denominator;
676 else
677 return 0;
678}
679
680double AliJCorrelations::GetGeoAccCorrIncl(double deltaEta){
681 // histo filler
682
683 if(fhistos->fhDEtaNearMixFromFile[fCentralityBin][fpttBin][fptaBin]->GetEntries()<1000) return GetGeoAccCorrFlat(deltaEta);
684
685 int bin = fhistos->fhDEtaNearMixFromFile[fCentralityBin][fpttBin][fptaBin]->FindBin(deltaEta);
686 double denominator = fhistos->fhDEtaNearMixFromFile[fCentralityBin][fpttBin][fptaBin]->GetBinContent(bin);
687 if(denominator > 1e-6)
688 return 1.0/denominator;
689 else
690 return 0;
691
692}
693
694double AliJCorrelations::DeltaPhi(double phi1, double phi2) {
695 // dphi
696 double res = atan2(sin(phi1-phi2), cos(phi1-phi2));
697 //double res = phi1 - phi2;
698 //return res>-kJPi/3.0 ? res : kJTwoPi+res ;
699 return res>-kJPi*9./20. ? res : kJTwoPi+res ;
700}
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721