PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHMEC.cxx

   1 // $Id$
   2
   3 //////////
   4 //Measure Jet-hadron correlations
   5 //Does event Mixing using AliEventPoolManager
   6 /////////
   7
   8 #include "AliAnalysisTaskEmcalJetHMEC.h"
   9
  10 #include "TChain.h"
  11 #include "TTree.h"
  12 #include "TList.h"
  13 #include "TH1F.h"
  14 #include "TH2F.h"
  15 #include "THnSparse.h"
  16 #include "TCanvas.h"
  17 #include <TClonesArray.h>
  18 #include <TParticle.h>
  19 #include "AliVTrack.h"
  20 #include "TParameter.h"
  21
  22 #include "AliAODEvent.h"
  23 #include "AliAnalysisTask.h"
  24 #include "AliAnalysisManager.h"
  25
  26 #include "AliESDEvent.h"
  27 #include "AliESDInputHandler.h"
  28 #include "AliESDCaloCluster.h"
  29 #include "AliESDVertex.h"
  30 #include "AliCentrality.h"
  31 #include "AliAODJet.h"
  32 #include "AliEmcalJet.h"
  33 #include "AliESDtrackCuts.h"
  34
  35 #include "TVector3.h"
  36 #include "AliPicoTrack.h"
  37 #include "AliEventPoolManager.h"
  38
  39 ClassImp(AliAnalysisTaskEmcalJetHMEC)
  40
  41 //________________________________________________________________________
  42 AliAnalysisTaskEmcalJetHMEC::AliAnalysisTaskEmcalJetHMEC() :
  43   AliAnalysisTaskEmcalJet("HMEC",kFALSE),
  44   fTracksName(""),
  45   fJetsName(""),
  46   fPhimin(-10),
  47   fPhimax(10),
  48   fEtamin(-0.9),
  49   fEtamax(0.9),
  50   fAreacut(0.0),
  51   fTrkBias(5),
  52   fClusBias(5),
  53   fTrkEta(0.9),
  54   fDoEventMixing(0),
  55   fMixingTracks(50000), fNMIXtracks(5000), fNMIXevents(5),
  56   fTriggerEventType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
  57   fDoEffCorrection(0), fEffFunctionCorrection(0),
  58   fDoLessSparseAxes(0), fDoWiderTrackBin(0),
  59   fCentBinSize(1),
  60   fESD(0),
  61   fAOD(0),
  62   fPoolMgr(0x0),
  63   fHistTrackPt(0),
  64   fHistCentrality(0),
  65   fHistJetEtaPhi(0),
  66   fHistJetHEtaPhi(0),
  67   fhnMixedEvents(0x0),
  68   fhnJH(0x0)
  69 {
  70   // Default Constructor
  71
  72   for(Int_t ipta=0; ipta<7; ipta++){
  73     fHistTrackEtaPhi[ipta]=0;
  74   }
  75
  76   for(Int_t icent = 0; icent<6; ++icent){
  77     fHistJetPt[icent]=0;
  78     fHistJetPtBias[icent]=0;
  79     fHistLeadJetPt[icent]=0;
  80     fHistLeadJetPtBias[icent]=0;
  81     fHistJetPtTT[icent]=0;
  82     for(Int_t iptjet = 0; iptjet<5; ++iptjet){
  83       for(Int_t ieta = 0; ieta<3; ++ieta){
  84             fHistJetH[icent][iptjet][ieta]=0;
  85             fHistJetHBias[icent][iptjet][ieta]=0;
  86             fHistJetHTT[icent][iptjet][ieta]=0;
  87       }
  88     }
  89   }
  90
  91 }
  92 //________________________________________________________________________
  93 AliAnalysisTaskEmcalJetHMEC::AliAnalysisTaskEmcalJetHMEC(const char *name) :
  94   AliAnalysisTaskEmcalJet(name,kTRUE),
  95   fTracksName(""),
  96   fJetsName(""),
  97   fPhimin(-10),
  98   fPhimax(10),
  99   fEtamin(-0.9),
 100   fEtamax(0.9),
 101   fAreacut(0.0),
 102   fTrkBias(5),
 103   fClusBias(5),
 104   fTrkEta(0.9),
 105   fDoEventMixing(0),
 106   fMixingTracks(50000), fNMIXtracks(5000), fNMIXevents(5),
 107   fTriggerEventType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
 108   fDoEffCorrection(0), fEffFunctionCorrection(0),
 109   fDoLessSparseAxes(0), fDoWiderTrackBin(0),
 110   fCentBinSize(1),
 111   fESD(0),
 112   fAOD(0),
 113   fPoolMgr(0x0),
 114   fHistTrackPt(0),
 115   fHistCentrality(0),
 116   fHistJetEtaPhi(0),
 117   fHistJetHEtaPhi(0),
 118   fhnMixedEvents(0x0),
 119   fhnJH(0x0)
 120 {
 121   // Constructor
 122   for(Int_t ipta=0; ipta<7; ipta++){
 123     fHistTrackEtaPhi[ipta]=0;
 124   }
 125   for(Int_t icent = 0; icent<6; ++icent){
 126     fHistJetPt[icent]=0;
 127     fHistJetPtBias[icent]=0;
 128     fHistLeadJetPt[icent]=0;
 129     fHistLeadJetPtBias[icent]=0;
 130     fHistJetPtTT[icent]=0;
 131     for(Int_t iptjet = 0; iptjet<5; ++iptjet){
 132       for(Int_t ieta = 0; ieta<3; ++ieta){
 133             fHistJetH[icent][iptjet][ieta]=0;
 134             fHistJetHBias[icent][iptjet][ieta]=0;
 135             fHistJetHTT[icent][iptjet][ieta]=0;
 136       }
 137     }
 138   }
 139
 140 }
 141
 142 //________________________________________________________________________
 143 void AliAnalysisTaskEmcalJetHMEC::UserCreateOutputObjects() {
 144   // Called once
 145   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 146   OpenFile(1);
 147
 148   // Create histograms
 149   fHistTrackPt = new TH1F("fHistTrackPt", "P_{T} distribution", 1000, 0.0, 100.0);
 150   fHistCentrality = new TH1F("fHistCentrality","centrality",100,0,100);
 151   fHistJetEtaPhi = new TH2F("fHistJetEtaPhi","Jet eta-phi",900,-1.8,1.8,720,-3.2,3.2);
 152   fHistJetHEtaPhi = new TH2F("fHistJetHEtaPhi","Jet-Hadron deta-dphi",900,-1.8,1.8,720,-1.6,4.8);
 153
 154   TString name;
 155
 156   for(Int_t ipta=0; ipta<7; ++ipta){
 157     name = Form("fHistTrackEtaPhi_%i", ipta);
 158     fHistTrackEtaPhi[ipta] = new TH2F(name,name,400,-1,1,720,0.0,2.0*TMath::Pi());
 159     fOutput->Add(fHistTrackEtaPhi[ipta]);
 160   }
 161
 162   for(Int_t icent = 0; icent<6; ++icent){
 163     name = Form("fHistJetPt_%i",icent);
 164     fHistJetPt[icent] = new TH1F(name,name,200,0,200);
 165     fOutput->Add(fHistJetPt[icent]);
 166
 167     name = Form("fHistJetPtBias_%i",icent);
 168     fHistJetPtBias[icent] = new TH1F(name,name,200,0,200);
 169     fOutput->Add(fHistJetPtBias[icent]);
 170
 171     name = Form("fHistLeadJetPt_%i",icent);
 172     fHistLeadJetPt[icent] = new TH1F(name,name,200,0,200);
 173     fOutput->Add(fHistLeadJetPt[icent]);
 174
 175     name = Form("fHistLeadJetPtBias_%i",icent);
 176     fHistLeadJetPtBias[icent] = new TH1F(name,name,200,0,200);
 177     fOutput->Add(fHistLeadJetPtBias[icent]);
 178
 179     name = Form("fHistJetPtTT_%i",icent);
 180     fHistJetPtTT[icent] = new TH1F(name,name,200,0,200);
 181     fOutput->Add(fHistJetPtTT[icent]);
 182
 183     for(Int_t iptjet = 0; iptjet<5; ++iptjet){
 184       for(Int_t ieta = 0; ieta<3; ++ieta){
 185             name = Form("fHistJetH_%i_%i_%i",icent,iptjet,ieta);
 186             fHistJetH[icent][iptjet][ieta]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
 187             fOutput->Add(fHistJetH[icent][iptjet][ieta]);
 188
 189         name = Form("fHistJetHBias_%i_%i_%i",icent,iptjet,ieta);
 190             fHistJetHBias[icent][iptjet][ieta]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
 191             fOutput->Add(fHistJetHBias[icent][iptjet][ieta]);
 192
 193             name = Form("fHistJetHTT_%i_%i_%i",icent,iptjet,ieta);
 194             fHistJetHTT[icent][iptjet][ieta]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
 195             fOutput->Add(fHistJetHTT[icent][iptjet][ieta]);
 196
 197       }
 198     }
 199   }
 200
 201   UInt_t cifras = 0; // bit coded, see GetDimParams() below
 202   if(fDoLessSparseAxes) {
 203     cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
 204   } else {
 205     cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
 206   }
 207   fhnJH = NewTHnSparseF("fhnJH", cifras);
 208   fhnJH->Sumw2();
 209   fOutput->Add(fhnJH);
 210
 211   if(fDoEventMixing){
 212     if(fDoLessSparseAxes) {
 213       cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
 214     } else {
 215       cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
 216     }
 217     fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
 218     fhnMixedEvents->Sumw2();
 219     fOutput->Add(fhnMixedEvents);
 220   }
 221
 222   fOutput->Add(fHistTrackPt);
 223   fOutput->Add(fHistCentrality);
 224   fOutput->Add(fHistJetEtaPhi);
 225   fOutput->Add(fHistJetHEtaPhi);
 226
 227   PostData(1, fOutput);
 228
 229   //Event Mixing
 230   Int_t trackDepth = fMixingTracks;
 231   Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
 232
 233   Int_t nZvtxBins  = 5+1+5;
 234   // bins for second buffer are shifted by 100 cm
 235   Double_t vertexBins[] = { -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, };
 236   Double_t* zvtxbin = vertexBins;
 237
 238 //  Int_t nCentralityBins  = 100;
 239   Int_t nCentralityBins = 100;
 240   Double_t mult = 1.0;
 241   if(fCentBinSize==1) {
 242     nCentralityBins = 100;
 243     mult = 1.0;
 244   } else if(fCentBinSize==2){
 245     nCentralityBins = 50;
 246     mult = 2.0;
 247   } else if(fCentBinSize==5){
 248     nCentralityBins = 20;
 249     mult = 5.0;
 250   } else if(fCentBinSize==10){
 251     nCentralityBins = 10;
 252     mult = 10.0;
 253   }
 254   Double_t centralityBins[nCentralityBins];
 255   for(Int_t ic=0; ic<nCentralityBins; ic++){
 256    centralityBins[ic]=mult*ic;
 257   }
 258
 259   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
 260
 261 }
 262
 263 //________________________________________________________________________
 264
 265 Double_t AliAnalysisTaskEmcalJetHMEC:: RelativePhi(Double_t mphi,Double_t vphi) {
 266
 267   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
 268   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
 269   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
 270   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
 271   Double_t dphi = vphi-mphi;
 272   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
 273   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
 274
 275   return dphi;//dphi in [-Pi, Pi]
 276 }
 277
 278 //________________________________________________________________________
 279 Int_t AliAnalysisTaskEmcalJetHMEC::GetCentBin(Double_t cent) const {
 280   // Get centrality bin.
 281
 282   Int_t centbin = -1;
 283   if (cent>=0 && cent<10) centbin = 0;
 284   else if (cent>=10 && cent<20) centbin = 1;
 285   else if (cent>=20 && cent<30) centbin = 2;
 286   else if (cent>=30 && cent<40) centbin = 3;
 287   else if (cent>=40 && cent<50) centbin = 4;
 288   else if (cent>=50 && cent<90) centbin = 5;
 289   return centbin;
 290 }
 291
 292 //________________________________________________________________________
 293 Int_t AliAnalysisTaskEmcalJetHMEC::GetEtaBin(Double_t eta) const {
 294   // Get eta bin for histos.
 295
 296   Int_t etabin = -1;
 297   if (TMath::Abs(eta)<=0.4) etabin = 0;
 298   else if (TMath::Abs(eta)>0.4 && TMath::Abs(eta)<0.8) etabin = 1;
 299   else if (TMath::Abs(eta)>=0.8) etabin = 2;
 300   return etabin;
 301 }
 302 //________________________________________________________________________
 303 Int_t AliAnalysisTaskEmcalJetHMEC::GetpTjetBin(Double_t pt) const
 304 {
 305   // Get jet pt  bin for histos.
 306
 307   Int_t ptbin = -1;
 308   if (pt>=15 && pt<20) ptbin = 0;
 309   else if (pt>=20 && pt<25) ptbin = 1;
 310   else if (pt>=25 && pt<30) ptbin = 2;
 311   else if (pt>=30 && pt<60) ptbin = 3;
 312   else if (pt>=60) ptbin = 4;
 313
 314   return ptbin;
 315 }
 316
 317 //________________________________________________________________________
 318 void AliAnalysisTaskEmcalJetHMEC::ExecOnce() {
 319   AliAnalysisTaskEmcalJet::ExecOnce();
 320
 321 }
 322
 323 //________________________________________________________________________
 324 Bool_t AliAnalysisTaskEmcalJetHMEC::Run() {
 325  // Main loop called for each event
 326   if(!fTracks){
 327     AliError(Form("No fTracks object!!\n"));
 328     return kTRUE;
 329   }
 330   if(!fJets){
 331     AliError(Form("No fJets object!!\n"));
 332     return kTRUE;
 333   }
 334
 335   // what kind of event do we have: AOD or ESD?
 336   Bool_t esdMode = kTRUE;
 337   if (dynamic_cast<AliAODEvent*>(InputEvent())) esdMode = kFALSE;
 338
 339   // if we have ESD event, set up ESD object
 340   if(esdMode){
 341     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 342     if (!fESD) {
 343       AliError(Form("ERROR: fESD not available\n"));
 344       return kTRUE;
 345     }
 346   }
 347
 348   // if we have AOD event, set up AOD object
 349   if(!esdMode){
 350     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 351     if(!fAOD) {
 352       AliError(Form("ERROR: fAOD not available\n"));
 353       return kTRUE;
 354     }
 355   }
 356
 357   TList *list = InputEvent()->GetList();
 358   if(!list) {
 359     AliError(Form("ERROR: list not attached\n"));
 360     return kTRUE;
 361   }
 362
 363   // get centrality
 364   if (fCent<0) {
 365     AliError(Form("Centrality negative: %f", fCent));
 366     return kTRUE;
 367   }
 368
 369   Int_t centbin = GetCentBin(fCent);
 370   if(centbin<0) return kTRUE;
 371
 372   Double_t fvertex[3]={0,0,0};
 373   InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
 374   Double_t zVtx=fvertex[2];
 375
 376   if(fabs(zVtx)>10.0) return kTRUE;
 377
 378   fHistCentrality->Fill(fCent);
 379
 380   TClonesArray *jets = 0;
 381   TClonesArray *tracks = 0;
 382
 383   tracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracks));
 384   if (!tracks) {
 385     AliError(Form("Pointer to tracks %s == 0", fTracks->GetName() ));
 386     return kTRUE;
 387   }
 388   const Int_t Ntracks=tracks->GetEntries();
 389
 390   jets= dynamic_cast<TClonesArray*>(list->FindObject(fJets));
 391   if (!jets) {
 392     AliError(Form("Pointer to tracks %s == 0", fJets->GetName() ));
 393     return kTRUE;
 394   }
 395   const Int_t Njets = jets->GetEntries();
 396
 397   //Leticia's loop to find hardest track
 398   Int_t iTT=-1;
 399   Double_t ptmax=-10;
 400
 401   for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++) {
 402     AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 403     if (!track) {
 404       printf("ERROR: Could not receive track %d\n", iTracks);
 405       continue;
 406     }
 407
 408     if(TMath::Abs(track->Eta())>0.9) continue;
 409     if(track->Pt()<0.15) continue;
 410     //iCount++;
 411     if(track->Pt()>ptmax){
 412           ptmax=track->Pt();
 413           iTT=iTracks;
 414     }
 415   }
 416
 417   Int_t ijethi=-1;
 418   Double_t highestjetpt=0.0;
 419   Int_t passedTTcut=0;
 420
 421   for (Int_t ijets = 0; ijets < Njets; ijets++){
 422       AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijets));
 423       if (!jet) continue;
 424       if(!AcceptthisJet(jet)) continue;
 425
 426       Double_t jetPt = jet->Pt();
 427
 428       if(highestjetpt<jetPt){
 429             ijethi=ijets;
 430             highestjetpt=jetPt;
 431       }
 432   }
 433
 434   // see if event is selected
 435   UInt_t trig = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 436
 437   for (Int_t ijet = 0; ijet < Njets; ijet++){
 438     AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
 439     if (!jet) continue;
 440
 441     // see if event is selected and our jet came from trigger event
 442     if (!(trig & fTriggerEventType)) continue;
 443     if (jet->Pt()<0.1) continue;
 444
 445     if(!AcceptthisJet(jet)) continue;
 446
 447     Double_t jetphi = jet->Phi();
 448     Double_t jetPt = jet->Pt();
 449     Double_t jeteta=jet->Eta();
 450
 451     Double_t leadjet=0;
 452     if (ijet==ijethi) leadjet=1;
 453
 454     fHistJetPt[centbin]->Fill(jet->Pt());
 455     fHistLeadJetPt[centbin]->Fill(jet->Pt());
 456
 457     if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
 458       fHistJetPtBias[centbin]->Fill(jet->Pt());
 459       fHistLeadJetPtBias[centbin]->Fill(jet->Pt());
 460     }
 461
 462     fHistJetEtaPhi->Fill(jet->Eta(),jetphi);
 463
 464     if(iTT>0){
 465           AliVTrack* TT = static_cast<AliVTrack*>(tracks->At(iTT));
 466           if(TMath::Abs(jetphi-TT->Phi()-TMath::Pi())<0.6) passedTTcut=1;
 467           else passedTTcut=0;
 468     }
 469
 470     if(passedTTcut)
 471           fHistJetPtTT[centbin]->Fill(jet->Pt());
 472
 473       Int_t iptjet=-1;
 474       iptjet=GetpTjetBin(jetPt);
 475       if(iptjet<0) continue;
 476
 477     if (jetPt > 15) {
 478
 479       for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++) {
 480             AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 481             if (!track) {
 482               printf("ERROR: Could not receive track %d\n", iTracks);
 483               continue;
 484             }
 485
 486             if(TMath::Abs(track->Eta())>fTrkEta) continue;
 487
 488             fHistTrackPt->Fill(track->Pt());
 489
 490             if (track->Pt()<0.15) continue;
 491
 492             Double_t trackphi = track->Phi();
 493             if (trackphi > TMath::Pi())
 494               trackphi = trackphi-2*TMath::Pi();
 495
 496             Double_t tracketa=track->Eta();
 497             Double_t trackpt=track->Pt();
 498             Double_t deta=tracketa-jeteta;
 499             Int_t ieta=GetEtaBin(deta);
 500             if (ieta<0) {
 501               AliError(Form("Eta Bin negative: %f", deta));
 502               continue;
 503             }
 504
 505             //Jet pt, track pt, dPhi,deta,fCent
 506             Double_t dphijh = RelativePhi(jetphi,trackphi);
 507             if (dphijh < -0.5*TMath::Pi())
 508               dphijh+= 2*TMath::Pi();
 509             if (dphijh > 1.5*TMath::Pi())
 510           dphijh-=2.*TMath::Pi();
 511
 512             fHistJetH[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
 513             fHistJetHEtaPhi->Fill(deta,dphijh);
 514
 515         // calculate single particle tracking efficiency for correlations
 516         Double_t trefficiency = -999;
 517         trefficiency = EffCorrection(track->Eta(), track->Pt(), fDoEffCorrection);
 518
 519             Double_t dR=sqrt(deta*deta+dphijh*dphijh);
 520
 521             if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
 522               fHistJetHBias[centbin][iptjet][ieta]->Fill(dphijh,trackpt);
 523
 524           if(fDoLessSparseAxes) { // check if we want all dimensions
 525             Double_t triggerEntries[6] = {fCent,jetPt,trackpt,deta,dphijh,leadjet};
 526             fhnJH->Fill(triggerEntries, 1.0/trefficiency);
 527           } else {
 528                 Double_t triggerEntries[8] = {fCent,jetPt,trackpt,deta,dphijh,leadjet,0.0,dR};
 529             fhnJH->Fill(triggerEntries, 1.0/trefficiency);
 530           }
 531             }
 532
 533             if(passedTTcut)
 534               fHistJetHTT[centbin][iptjet][ieta]->Fill(dphijh,trackpt);
 535
 536           } //track loop
 537     }//jet pt cut
 538   }//jet loop
 539
 540   //Prepare to do event mixing
 541
 542   // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
 543   TObjArray* tracksClone = 0x0;
 544
 545   if(fDoEventMixing>0){
 546
 547     // event mixing
 548
 549     // 1. First get an event pool corresponding in mult (cent) and
 550     //    zvertex to the current event. Once initialized, the pool
 551     //    should contain nMix (reduced) events. This routine does not
 552     //    pre-scan the chain. The first several events of every chain
 553     //    will be skipped until the needed pools are filled to the
 554     //    specified depth. If the pool categories are not too rare, this
 555     //    should not be a problem. If they are rare, you could lose
 556     //    statistics.
 557
 558     // 2. Collect the whole pool's content of tracks into one TObjArray
 559     //    (bgTracks), which is effectively a single background super-event.
 560
 561     // 3. The reduced and bgTracks arrays must both be passed into
 562     //    FillCorrelations(). Also nMix should be passed in, so a weight
 563     //    of 1./nMix can be applied.
 564
 565     UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 566     // if event was not selected (triggered) for any reseason (should never happen) then return
 567     if (trigger==0)  return kTRUE;
 568
 569     AliEventPool* pool = fPoolMgr->GetEventPool(fCent, zVtx);
 570
 571     if (!pool){
 572       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCent, zVtx));
 573       return kTRUE;
 574     }
 575
 576 if(trigger & fTriggerEventType) {
 577     //check for a trigger jet
 578     if (pool->IsReady() || pool->NTracksInPool() > fNMIXtracks || pool->GetCurrentNEvents() >= fNMIXevents) {
 579
 580         for (Int_t ijet = 0; ijet < Njets; ijet++){
 581           Double_t leadjet=0;
 582           if (ijet==ijethi) leadjet=1;
 583
 584           AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
 585       if (!jet) continue;
 586
 587           if(!AcceptthisJet(jet)) continue;
 588
 589           Double_t jetPt = jet->Pt();
 590           Double_t jetphi = jet->Phi();
 591           Double_t jeteta=jet->Eta();
 592
 593       // make sure event contains jet above our threshold (reduce stats of sparse)
 594       if (jetPt < 15) continue;
 595
 596           Int_t nMix = pool->GetCurrentNEvents();
 597
 598           //Fill for biased jet triggers only
 599           if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
 600
 601             // Fill mixed-event histos here
 602             for (Int_t jMix=0; jMix<nMix; jMix++) {
 603                   TObjArray* bgTracks = pool->GetEvent(jMix);
 604                   const Int_t Nbgtrks = bgTracks->GetEntries();
 605
 606           for(Int_t ibg=0; ibg<Nbgtrks; ibg++){
 607                     AliPicoTrack *part = static_cast<AliPicoTrack*>(bgTracks->At(ibg));
 608                     if(!part) continue;
 609             if(TMath::Abs(part->Eta())>0.9) continue;
 610             if(part->Pt()<0.15) continue;
 611
 612                     Double_t DEta = part->Eta()-jeteta;
 613                     Double_t DPhi = RelativePhi(jetphi,part->Phi());
 614
 615             // calculate single particle tracking efficiency of mixed events for correlations
 616             Double_t mixefficiency = -999;
 617             mixefficiency = EffCorrection(part->Eta(), part->Pt(), fDoEffCorrection);
 618
 619                     Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
 620                     if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
 621                     if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
 622             if(fDoLessSparseAxes) {  // check if we want all the axis filled
 623               Double_t triggerEntries[6] = {fCent,jetPt,part->Pt(),DEta,DPhi,leadjet};
 624               fhnMixedEvents->Fill(triggerEntries,1./(nMix*mixefficiency));
 625             } else {
 626                       Double_t triggerEntries[8] = {fCent,jetPt,part->Pt(),DEta,DPhi,leadjet,0.0,DR};
 627               fhnMixedEvents->Fill(triggerEntries,1./(nMix*mixefficiency));
 628                     }
 629                     }
 630               }
 631             }
 632           }
 633     }
 634 }
 635
 636     if(trigger & fMixingEventType) {
 637       tracksClone = CloneAndReduceTrackList(tracks);
 638
 639       //update pool if jet in event or not
 640       pool->UpdatePool(tracksClone);
 641     }
 642
 643   } // end of event mixing
 644
 645   return kTRUE;
 646 }
 647
 648 //________________________________________________________________________
 649 void AliAnalysisTaskEmcalJetHMEC::Terminate(Option_t *)
 650 {
 651   //just terminate
 652 }
 653
 654 //________________________________________________________________________
 655 Int_t AliAnalysisTaskEmcalJetHMEC::AcceptthisJet(AliEmcalJet *jet)
 656 {
 657   //applies all jet cuts except pt
 658   float jetphi = jet->Phi();
 659   if (jetphi>TMath::Pi())
 660     jetphi = jetphi-2*TMath::Pi();
 661
 662   if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax))
 663     return 0;
 664   if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax))
 665     return 0;
 666   if (jet->Area()<fAreacut)
 667     return 0;
 668   //prevents 0 area jets from sneaking by when area cut == 0
 669   if (jet->Area()==0)
 670     return 0;
 671   //exclude jets with extremely high pt tracks which are likely misreconstructed
 672   if(jet->MaxTrackPt()>100)
 673     return 0;
 674
 675   //passed all above cuts
 676   return 1;
 677 }
 678
 679 //________________________________________________________________________
 680
 681 THnSparse* AliAnalysisTaskEmcalJetHMEC::NewTHnSparseF(const char* name, UInt_t entries){
 682    // generate new THnSparseF, axes are defined in GetDimParams()
 683
 684    Int_t count = 0;
 685    UInt_t tmp = entries;
 686    while(tmp!=0){
 687       count++;
 688       tmp = tmp &~ -tmp;  // clear lowest bit
 689    }
 690
 691    TString hnTitle(name);
 692    const Int_t dim = count;
 693    Int_t nbins[dim];
 694    Double_t xmin[dim];
 695    Double_t xmax[dim];
 696
 697    Int_t i=0;
 698    Int_t c=0;
 699    while(c<dim && i<32){
 700       if(entries&(1<<i)){
 701
 702          TString label("");
 703          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
 704          hnTitle += Form(";%s",label.Data());
 705          c++;
 706       }
 707
 708       i++;
 709    }
 710    hnTitle += ";";
 711
 712    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
 713 }
 714
 715 void AliAnalysisTaskEmcalJetHMEC::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
 716 {
 717    // stores label and binning of axis for THnSparse
 718
 719    const Double_t pi = TMath::Pi();
 720
 721    switch(iEntry){
 722
 723    case 0:
 724       label = "V0 centrality (%)";
 725          nbins = 10;
 726          xmin = 0.;
 727          xmax = 100.;
 728          break;
 729
 730    case 1:
 731       label = "corrected jet pt";
 732       nbins = 20;
 733       xmin = 0.;
 734       xmax = 200.;
 735       break;
 736
 737    case 2:
 738       if(fDoWiderTrackBin) {
 739         label = "track pT";
 740         nbins = 40;
 741         xmin = 0.;
 742         xmax = 10.;
 743       } else {
 744         label = "track pT";
 745         nbins = 100;
 746         xmin = 0.;
 747         xmax = 10;
 748       }
 749       break;
 750
 751    case 3:
 752       label = "deltaEta";
 753       nbins = 24;
 754       xmin = -1.2;
 755       xmax = 1.2;
 756       break;
 757
 758    case 4:
 759       label = "deltaPhi";
 760       nbins = 72;
 761       xmin = -0.5*pi;
 762       xmax = 1.5*pi;
 763       break;
 764
 765     case 5:
 766       label = "leading jet";
 767       nbins = 3;
 768       xmin = -0.5;
 769       xmax = 2.5;
 770       break;
 771
 772     case 6:
 773       label = "trigger track";
 774       nbins =10;
 775       xmin = 0;
 776       xmax = 50;
 777       break;
 778
 779     case 7:
 780       label = "deltaR";
 781       nbins = 10;
 782       xmin = 0.;
 783       xmax = 5.0;
 784       break;
 785
 786     case 8:
 787       label = "leading track";
 788       nbins = 13;
 789       xmin = 0;
 790       xmax = 50;
 791       break;
 792    }
 793 }
 794
 795 //_________________________________________________
 796 // From CF event mixing code PhiCorrelations
 797 TObjArray* AliAnalysisTaskEmcalJetHMEC::CloneAndReduceTrackList(TObjArray* tracks)
 798 {
 799   // clones a track list by using AliPicoTrack which uses much less memory (used for event mixing)
 800
 801   TObjArray* tracksClone = new TObjArray;
 802   tracksClone->SetOwner(kTRUE);
 803
 804   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
 805   {
 806     AliVParticle* particle = (AliVParticle*) tracks->At(i);
 807     if(TMath::Abs(particle->Eta())>fTrkEta) continue;
 808     if(particle->Pt()<0.15) continue;
 809
 810     Double_t trackpt=particle->Pt();
 811
 812     Int_t hadbin=-1;
 813     if(trackpt<0.5) hadbin=0;
 814     else if(trackpt<1) hadbin=1;
 815     else if(trackpt<2) hadbin=2;
 816     else if(trackpt<3) hadbin=3;
 817     else if(trackpt<5) hadbin=4;
 818     else if(trackpt<8) hadbin=5;
 819     else if(trackpt<20) hadbin=6;
 820
 821     if(hadbin>-1) fHistTrackEtaPhi[hadbin]->Fill(particle->Eta(),particle->Phi());
 822
 823     tracksClone->Add(new AliPicoTrack(particle->Pt(), particle->Eta(), particle->Phi(), particle->Charge(), 0, 0, 0, 0));
 824   }
 825
 826   return tracksClone;
 827 }
 828
 829 //________________________________________________________________________
 830 Double_t AliAnalysisTaskEmcalJetHMEC::EffCorrection(Double_t trackETA, Double_t trackPT, Int_t effSwitch) const {
 831   // x-variable = track pt, y-variable = track eta
 832   Double_t x = trackPT;
 833   Double_t y = trackETA;
 834   Double_t TRefficiency = -999;
 835   Int_t runNUM = fCurrentRunNumber;
 836   Int_t runSwitchGood = -999;
 837   Int_t centbin = -99;
 838   Double_t etaaxis = 0;
 839   Double_t ptaxis = 0;
 840
 841   if(effSwitch < 1) {
 842     if ((runNUM == 169975 || runNUM == 169981 || runNUM == 170038 || runNUM == 170040 || runNUM == 170083 || runNUM == 170084 || runNUM == 170085 || runNUM == 170088 || runNUM == 170089 || runNUM == 170091 || runNUM == 170155 || runNUM == 170159 || runNUM == 170163 || runNUM == 170193 || runNUM == 170195 || runNUM == 170203 || runNUM == 170204 || runNUM == 170228 || runNUM == 170230 || runNUM == 170268 || runNUM == 170269 || runNUM == 170270 || runNUM == 170306 || runNUM == 170308 || runNUM == 170309)) runSwitchGood = 0;
 843
 844     if ((runNUM == 167903 || runNUM == 167915 || runNUM == 167987 || runNUM == 167988 || runNUM == 168066 || runNUM == 168068 || runNUM == 168069 || runNUM == 168076 || runNUM == 168104 || runNUM == 168107 || runNUM == 168108 || runNUM == 168115 || runNUM == 168212 || runNUM == 168310 || runNUM == 168311 || runNUM == 168322 || runNUM == 168325 || runNUM == 168341 || runNUM == 168342 || runNUM == 168361 || runNUM == 168362 || runNUM == 168458 || runNUM == 168460 || runNUM == 168461 || runNUM == 168464 || runNUM == 168467 || runNUM == 168511 || runNUM == 168512 || runNUM == 168777 || runNUM == 168826 || runNUM == 168984 || runNUM == 168988 || runNUM == 168992 || runNUM == 169035 || runNUM == 169091 || runNUM == 169094 || runNUM == 169138 || runNUM == 169143 || runNUM == 169144 || runNUM == 169145 || runNUM == 169148 || runNUM == 169156 || runNUM == 169160 || runNUM == 169167 || runNUM == 169238 || runNUM == 169411 || runNUM == 169415 || runNUM == 169417 || runNUM == 169835 || runNUM == 169837 || runNUM == 169838 || runNUM == 169846 || runNUM == 169855 || runNUM == 169858 || runNUM == 169859 || runNUM == 169923 || runNUM == 169956 || runNUM == 170027 || runNUM == 170036 || runNUM == 170081)) runSwitchGood = 1;
 845
 846     if(fCent>=0 && fCent<10) centbin = 0;
 847     else if (fCent>=10 && fCent<30)     centbin = 1;
 848     else if (fCent>=30 && fCent<50) centbin = 2;
 849     else if (fCent>=50 && fCent<90)     centbin = 3;
 850
 851     if(runSwitchGood == 0 && centbin == 0) effSwitch = 2;
 852     if(runSwitchGood == 0 && centbin == 1) effSwitch = 3;
 853     if(runSwitchGood == 0 && centbin == 2) effSwitch = 4;
 854     if(runSwitchGood == 0 && centbin == 3) effSwitch = 5;
 855     if(runSwitchGood == 1 && centbin == 0) effSwitch = 6;
 856     if(runSwitchGood == 1 && centbin == 1) effSwitch = 7;
 857     if(runSwitchGood == 1 && centbin == 2) effSwitch = 8;
 858     if(runSwitchGood == 1 && centbin == 3) effSwitch = 9;
 859
 860   }
 861
 862   // 0-10% centrality: Semi-Good Runs
 863   Double_t p0_10SG[17] = {9.15846e-01, 7.54031e-02, 1.11627e+00, -2.35529e-02, 8.03415e-01, 9.43327e-03, -3.29811e-04, 1.66757e+00, 4.64069e-02, 1.84170e-01, -2.09091e-01, 9.56263e-01, 3.55602e-01, -5.75471e-01, 1.12262e+00, 1.96927e-03, 4.23474e-01};
 864   // 10-30% centrality: Semi-Good Runs
 865   Double_t p10_30SG[17] = {9.27764e-01, 7.69155e-02, 1.11896e+00, -2.55034e-02, 7.57674e-01, 3.68652e-02, -3.75523e-03, 1.64452e+00, 4.49358e-02, 1.92339e-01, -1.99106e-01, 9.59782e-01, 3.50657e-01, -6.07217e-01, 1.09027e+00, 1.58922e-03, 4.57454e-01};
 866   // 30-50% centrality: Semi-Good Runs
 867   Double_t p30_50SG[17] = {9.99833e-01, 8.31729e-02, 1.21066e+00, -2.99180e-02, 8.89407e-01, 2.72632e-02, -6.47995e-03, 1.07181e+00, 3.52573e-03, 3.02167e-01, -8.80560e-02, 9.08866e-01, 3.02663e-01, -5.85203e-01, 9.99932e-01, 2.11319e-03, 5.17128e-01};
 868   // 50-90% centrality: Semi-Good Runs
 869   Double_t p50_90SG[17] = {9.99962e-01, 8.44121e-02, 1.24631e+00, -2.54355e-02, 7.28692e-01, 8.89584e-02, -1.02253e-02, 1.10620e+00, 3.80133e-03, 2.74054e-01, -9.15362e-02, 9.04696e-01, 2.84234e-01, -5.67995e-01, 9.99870e-01, 1.89310e-03, 4.91358e-01};
 870
 871   // 0-10% centrality: Good Runs
 872   Double_t p0_10G[17] = {9.41769e-01, 7.67433e-02, 1.13335e+00, -2.66193e-02, 8.30263e-01, 5.19029e-03, 3.89549e-05, 1.64552e+00, 3.71978e-02, 1.92558e-01, -2.14009e-01, 9.49016e-01, 1.98726e-01, -2.77058e-01, 1.09062e+00, 1.90219e-03, 4.25480e-01};
 873   // 10-30% centrality: Good Runs
 874   Double_t p10_30G[17] = {9.60095e-01, 7.75887e-02, 1.12192e+00, -2.94920e-02, 8.28189e-01, 1.31810e-02, -1.20607e-03, 1.61751e+00, 3.60181e-02, 2.00906e-01, -1.99320e-01, 9.49605e-01, 1.77662e-01, -2.72506e-01, 1.05439e+00, 1.63010e-03, 4.66318e-01};
 875   // 30-50% centrality: Good Runs
 876   Double_t p30_50G[17] = {9.54634e-01, 8.21353e-02, 1.15674e+00, -3.30327e-02, 6.89883e-01, 8.50469e-02, -1.11555e-02, 1.23933e+00, 4.23219e-03, 2.80843e-01, -1.11685e-01, 9.88261e-01, 1.01834e-01, -1.91367e-01, 1.09119e+00, 1.63848e-03, 4.36947e-01};
 877   // 50-90% centrality: Good Runs
 878   Double_t p50_90G[17] = {9.70498e-01, 8.25935e-02, 1.15688e+00, -3.46965e-02, 7.25662e-01, 6.58188e-02, -7.40451e-03, 1.31635e+00, 5.30976e-03, 2.37158e-01, -1.21794e-01, 9.74182e-01, 1.08888e-01, -2.14606e-01, 1.07135e+00, 1.67715e-03, 4.47729e-01};
 879
 880   // set up a switch for different parameter values...
 881   switch(effSwitch) {
 882     case 1 :
 883       // first switch value - TRefficiency not used so = 1
 884       TRefficiency = 1.0;
 885       break;
 886
 887     case 2 :
 888       // Parameter values for Semi-GOOD TPC (LHC11h) runs (0-10%):
 889       ptaxis = (x<2.9)*(p0_10SG[0]*exp(-pow(p0_10SG[1]/x,p0_10SG[2])) + p0_10SG[3]*x) + (x>=2.9)*(p0_10SG[4] + p0_10SG[5]*x + p0_10SG[6]*x*x);
 890       etaaxis = (y<0.0)*(p0_10SG[7]*exp(-pow(p0_10SG[8]/TMath::Abs(y+0.9),p0_10SG[9])) + p0_10SG[10]*y) + (y>=0.0 && y<=0.4)*(p0_10SG[11] + p0_10SG[12]*y + p0_10SG[13]*y*y) + (y>0.4)*(p0_10SG[14]*exp(-pow(p0_10SG[15]/TMath::Abs(-y+0.9),p0_10SG[16])));
 891       TRefficiency = ptaxis*etaaxis;
 892       break;
 893
 894     case 3 :
 895       // Parameter values for Semi-GOOD TPC (LHC11h) runs (10-30%):
 896       ptaxis = (x<2.9)*(p10_30SG[0]*exp(-pow(p10_30SG[1]/x,p10_30SG[2])) + p10_30SG[3]*x) + (x>=2.9)*(p10_30SG[4] + p10_30SG[5]*x + p10_30SG[6]*x*x);
 897       etaaxis = (y<0.0)*(p10_30SG[7]*exp(-pow(p10_30SG[8]/TMath::Abs(y+0.9),p10_30SG[9])) + p10_30SG[10]*y) + (y>=0.0 && y<=0.4)*(p10_30SG[11] + p10_30SG[12]*y + p10_30SG[13]*y*y) + (y>0.4)*(p10_30SG[14]*exp(-pow(p10_30SG[15]/TMath::Abs(-y+0.9),p10_30SG[16])));
 898       TRefficiency = ptaxis*etaaxis;
 899       break;
 900
 901     case 4 :
 902       // Parameter values for Semi-GOOD TPC (LHC11h) runs (30-50%):
 903       ptaxis = (x<2.9)*(p30_50SG[0]*exp(-pow(p30_50SG[1]/x,p30_50SG[2])) + p30_50SG[3]*x) + (x>=2.9)*(p30_50SG[4] + p30_50SG[5]*x + p30_50SG[6]*x*x);
 904       etaaxis = (y<0.0)*(p30_50SG[7]*exp(-pow(p30_50SG[8]/TMath::Abs(y+0.9),p30_50SG[9])) + p30_50SG[10]*y) + (y>=0.0 && y<=0.4)*(p30_50SG[11] + p30_50SG[12]*y + p30_50SG[13]*y*y) + (y>0.4)*(p30_50SG[14]*exp(-pow(p30_50SG[15]/TMath::Abs(-y+0.9),p30_50SG[16])));
 905       TRefficiency = ptaxis*etaaxis;
 906       break;
 907
 908     case 5 :
 909       // Parameter values for Semi-GOOD TPC (LHC11h) runs (50-90%):
 910       ptaxis = (x<2.9)*(p50_90SG[0]*exp(-pow(p50_90SG[1]/x,p50_90SG[2])) + p50_90SG[3]*x) + (x>=2.9)*(p50_90SG[4] + p50_90SG[5]*x + p50_90SG[6]*x*x);
 911       etaaxis = (y<0.0)*(p50_90SG[7]*exp(-pow(p50_90SG[8]/TMath::Abs(y+0.9),p50_90SG[9])) + p50_90SG[10]*y) + (y>=0.0 && y<=0.4)*(p50_90SG[11] + p50_90SG[12]*y + p50_90SG[13]*y*y) + (y>0.4)*(p50_90SG[14]*exp(-pow(p50_90SG[15]/TMath::Abs(-y+0.9),p50_90SG[16])));
 912       TRefficiency = ptaxis*etaaxis;
 913       break;
 914
 915     case 6 :
 916       // Parameter values for GOOD TPC (LHC11h) runs (0-10%):
 917       ptaxis = (x<2.9)*(p0_10G[0]*exp(-pow(p0_10G[1]/x,p0_10G[2])) + p0_10G[3]*x) + (x>=2.9)*(p0_10G[4] + p0_10G[5]*x + p0_10G[6]*x*x);
 918       etaaxis = (y<0.0)*(p0_10G[7]*exp(-pow(p0_10G[8]/TMath::Abs(y+0.9),p0_10G[9])) + p0_10G[10]*y) + (y>=0.0 && y<=0.4)*(p0_10G[11] + p0_10G[12]*y + p0_10G[13]*y*y) + (y>0.4)*(p0_10G[14]*exp(-pow(p0_10G[15]/TMath::Abs(-y+0.9),p0_10G[16])));
 919       TRefficiency = ptaxis*etaaxis;
 920       break;
 921
 922     case 7 :
 923       // Parameter values for GOOD TPC (LHC11h) runs (10-30%):
 924       ptaxis = (x<2.9)*(p10_30G[0]*exp(-pow(p10_30G[1]/x,p10_30G[2])) + p10_30G[3]*x) + (x>=2.9)*(p10_30G[4] + p10_30G[5]*x + p10_30G[6]*x*x);
 925       etaaxis = (y<0.0)*(p10_30G[7]*exp(-pow(p10_30G[8]/TMath::Abs(y+0.9),p10_30G[9])) + p10_30G[10]*y) + (y>=0.0 && y<=0.4)*(p10_30G[11] + p10_30G[12]*y + p10_30G[13]*y*y) + (y>0.4)*(p10_30G[14]*exp(-pow(p10_30G[15]/TMath::Abs(-y+0.9),p10_30G[16])));
 926       TRefficiency = ptaxis*etaaxis;
 927       break;
 928
 929     case 8 :
 930       // Parameter values for GOOD TPC (LHC11h) runs (30-50%):
 931       ptaxis = (x<2.9)*(p30_50G[0]*exp(-pow(p30_50G[1]/x,p30_50G[2])) + p30_50G[3]*x) + (x>=2.9)*(p30_50G[4] + p30_50G[5]*x + p30_50G[6]*x*x);
 932       etaaxis = (y<0.0)*(p30_50G[7]*exp(-pow(p30_50G[8]/TMath::Abs(y+0.9),p30_50G[9])) + p30_50G[10]*y) + (y>=0.0 && y<=0.4)*(p30_50G[11] + p30_50G[12]*y + p30_50G[13]*y*y) + (y>0.4)*(p30_50G[14]*exp(-pow(p30_50G[15]/TMath::Abs(-y+0.9),p30_50G[16])));
 933       TRefficiency = ptaxis*etaaxis;
 934       break;
 935
 936     case 9 :
 937       // Parameter values for GOOD TPC (LHC11h) runs (50-90%):
 938       ptaxis = (x<2.9)*(p50_90G[0]*exp(-pow(p50_90G[1]/x,p50_90G[2])) + p50_90G[3]*x) + (x>=2.9)*(p50_90G[4] + p50_90G[5]*x + p50_90G[6]*x*x);
 939       etaaxis = (y<0.0)*(p50_90G[7]*exp(-pow(p50_90G[8]/TMath::Abs(y+0.9),p50_90G[9])) + p50_90G[10]*y) + (y>=0.0 && y<=0.4)*(p50_90G[11] + p50_90G[12]*y + p50_90G[13]*y*y) + (y>0.4)*(p50_90G[14]*exp(-pow(p50_90G[15]/TMath::Abs(-y+0.9),p50_90G[16])));
 940       TRefficiency = ptaxis*etaaxis;
 941       break;
 942
 943     default :
 944       // no Efficiency Switch option selected.. therefore don't correct, and set eff = 1
 945       TRefficiency = 1.0;
 946
 947   }
 948
 949   return TRefficiency;
 950 }
 951