update to master versions
[u/mrichter/AliRoot.git] / PWGCF / Correlations / JCORRAN / AliJCorrelations.cxx
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
37 AliJCorrelations::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
98 AliJCorrelations::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
135 AliJCorrelations::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
183 AliJCorrelations& 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
238 void 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 //=============================================================================================
246 void 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()==kJHadron && ftk2->GetParticleType()==kJHadron) 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
327 void 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
386 void 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
409 void 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
529 void 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
546 void 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
559 void 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
593 void 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
667 double 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
680 double 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
694 double 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