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()==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
 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