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   fDoLessSparseAxes(0), fDoWiderTrackBin(0),
  58   fCentBinSize(1),
  59   fESD(0),
  60   fAOD(0),
  61   fPoolMgr(0x0),
  62   fHistTrackPt(0),
  63   fHistCentrality(0),
  64   fHistJetEtaPhi(0),
  65   fHistJetHEtaPhi(0),
  66   fhnMixedEvents(0x0),
  67   fhnJH(0x0)
  68 {
  69   // Default Constructor
  70
  71   for(Int_t ipta=0; ipta<7; ipta++){
  72     fHistTrackEtaPhi[ipta]=0;
  73   }
  74
  75   for(Int_t icent = 0; icent<6; ++icent){
  76     fHistJetPt[icent]=0;
  77     fHistJetPtBias[icent]=0;
  78     fHistLeadJetPt[icent]=0;
  79     fHistLeadJetPtBias[icent]=0;
  80     fHistJetPtTT[icent]=0;
  81     for(Int_t iptjet = 0; iptjet<5; ++iptjet){
  82       for(Int_t ieta = 0; ieta<3; ++ieta){
  83             fHistJetH[icent][iptjet][ieta]=0;
  84             fHistJetHBias[icent][iptjet][ieta]=0;
  85             fHistJetHTT[icent][iptjet][ieta]=0;
  86       }
  87     }
  88   }
  89
  90 }
  91 //________________________________________________________________________
  92 AliAnalysisTaskEmcalJetHMEC::AliAnalysisTaskEmcalJetHMEC(const char *name) :
  93   AliAnalysisTaskEmcalJet(name,kTRUE),
  94   fTracksName(""),
  95   fJetsName(""),
  96   fPhimin(-10),
  97   fPhimax(10),
  98   fEtamin(-0.9),
  99   fEtamax(0.9),
 100   fAreacut(0.0),
 101   fTrkBias(5),
 102   fClusBias(5),
 103   fTrkEta(0.9),
 104   fDoEventMixing(0),
 105   fMixingTracks(50000), fNMIXtracks(5000), fNMIXevents(5),
 106   fTriggerEventType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
 107   fDoLessSparseAxes(0), fDoWiderTrackBin(0),
 108   fCentBinSize(1),
 109   fESD(0),
 110   fAOD(0),
 111   fPoolMgr(0x0),
 112   fHistTrackPt(0),
 113   fHistCentrality(0),
 114   fHistJetEtaPhi(0),
 115   fHistJetHEtaPhi(0),
 116   fhnMixedEvents(0x0),
 117   fhnJH(0x0)
 118 {
 119   // Constructor
 120   for(Int_t ipta=0; ipta<7; ipta++){
 121     fHistTrackEtaPhi[ipta]=0;
 122   }
 123   for(Int_t icent = 0; icent<6; ++icent){
 124     fHistJetPt[icent]=0;
 125     fHistJetPtBias[icent]=0;
 126     fHistLeadJetPt[icent]=0;
 127     fHistLeadJetPtBias[icent]=0;
 128     fHistJetPtTT[icent]=0;
 129     for(Int_t iptjet = 0; iptjet<5; ++iptjet){
 130       for(Int_t ieta = 0; ieta<3; ++ieta){
 131             fHistJetH[icent][iptjet][ieta]=0;
 132             fHistJetHBias[icent][iptjet][ieta]=0;
 133             fHistJetHTT[icent][iptjet][ieta]=0;
 134       }
 135     }
 136   }
 137
 138 }
 139
 140 //________________________________________________________________________
 141 void AliAnalysisTaskEmcalJetHMEC::UserCreateOutputObjects() {
 142   // Called once
 143   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 144   OpenFile(1);
 145
 146   // Create histograms
 147   fHistTrackPt = new TH1F("fHistTrackPt", "P_{T} distribution", 1000, 0.0, 100.0);
 148   fHistCentrality = new TH1F("fHistCentrality","centrality",100,0,100);
 149   fHistJetEtaPhi = new TH2F("fHistJetEtaPhi","Jet eta-phi",900,-1.8,1.8,720,-3.2,3.2);
 150   fHistJetHEtaPhi = new TH2F("fHistJetHEtaPhi","Jet-Hadron deta-dphi",900,-1.8,1.8,720,-1.6,4.8);
 151
 152   TString name;
 153
 154   for(Int_t ipta=0; ipta<7; ++ipta){
 155     name = Form("fHistTrackEtaPhi_%i", ipta);
 156     fHistTrackEtaPhi[ipta] = new TH2F(name,name,400,-1,1,720,0.0,2.0*TMath::Pi());
 157     fOutput->Add(fHistTrackEtaPhi[ipta]);
 158   }
 159
 160   for(Int_t icent = 0; icent<6; ++icent){
 161     name = Form("fHistJetPt_%i",icent);
 162     fHistJetPt[icent] = new TH1F(name,name,200,0,200);
 163     fOutput->Add(fHistJetPt[icent]);
 164
 165     name = Form("fHistJetPtBias_%i",icent);
 166     fHistJetPtBias[icent] = new TH1F(name,name,200,0,200);
 167     fOutput->Add(fHistJetPtBias[icent]);
 168
 169     name = Form("fHistLeadJetPt_%i",icent);
 170     fHistLeadJetPt[icent] = new TH1F(name,name,200,0,200);
 171     fOutput->Add(fHistLeadJetPt[icent]);
 172
 173     name = Form("fHistLeadJetPtBias_%i",icent);
 174     fHistLeadJetPtBias[icent] = new TH1F(name,name,200,0,200);
 175     fOutput->Add(fHistLeadJetPtBias[icent]);
 176
 177     name = Form("fHistJetPtTT_%i",icent);
 178     fHistJetPtTT[icent] = new TH1F(name,name,200,0,200);
 179     fOutput->Add(fHistJetPtTT[icent]);
 180
 181     for(Int_t iptjet = 0; iptjet<5; ++iptjet){
 182       for(Int_t ieta = 0; ieta<3; ++ieta){
 183             name = Form("fHistJetH_%i_%i_%i",icent,iptjet,ieta);
 184             fHistJetH[icent][iptjet][ieta]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
 185             fOutput->Add(fHistJetH[icent][iptjet][ieta]);
 186
 187         name = Form("fHistJetHBias_%i_%i_%i",icent,iptjet,ieta);
 188             fHistJetHBias[icent][iptjet][ieta]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
 189             fOutput->Add(fHistJetHBias[icent][iptjet][ieta]);
 190
 191             name = Form("fHistJetHTT_%i_%i_%i",icent,iptjet,ieta);
 192             fHistJetHTT[icent][iptjet][ieta]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
 193             fOutput->Add(fHistJetHTT[icent][iptjet][ieta]);
 194
 195       }
 196     }
 197   }
 198
 199   UInt_t cifras = 0; // bit coded, see GetDimParams() below
 200   if(fDoLessSparseAxes) {
 201     cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
 202   } else {
 203     cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
 204   }
 205   fhnJH = NewTHnSparseF("fhnJH", cifras);
 206   fhnJH->Sumw2();
 207   fOutput->Add(fhnJH);
 208
 209   if(fDoEventMixing){
 210     if(fDoLessSparseAxes) {
 211       cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
 212     } else {
 213       cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
 214     }
 215     fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
 216     fhnMixedEvents->Sumw2();
 217     fOutput->Add(fhnMixedEvents);
 218   }
 219
 220   fOutput->Add(fHistTrackPt);
 221   fOutput->Add(fHistCentrality);
 222   fOutput->Add(fHistJetEtaPhi);
 223   fOutput->Add(fHistJetHEtaPhi);
 224
 225   PostData(1, fOutput);
 226
 227   //Event Mixing
 228   Int_t trackDepth = fMixingTracks;
 229   Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
 230
 231   Int_t nZvtxBins  = 5+1+5;
 232   // bins for second buffer are shifted by 100 cm
 233   Double_t vertexBins[] = { -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, };
 234   Double_t* zvtxbin = vertexBins;
 235
 236 //  Int_t nCentralityBins  = 100;
 237   Int_t nCentralityBins = 100;
 238   Double_t mult = 1.0;
 239   if(fCentBinSize==1) {
 240     nCentralityBins = 100;
 241     mult = 1.0;
 242   } else if(fCentBinSize==2){
 243     nCentralityBins = 50;
 244     mult = 2.0;
 245   } else if(fCentBinSize==5){
 246     nCentralityBins = 20;
 247     mult = 5.0;
 248   } else if(fCentBinSize==10){
 249     nCentralityBins = 10;
 250     mult = 10.0;
 251   }
 252   Double_t centralityBins[nCentralityBins];
 253   for(Int_t ic=0; ic<nCentralityBins; ic++){
 254    centralityBins[ic]=mult*ic;
 255   }
 256
 257   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
 258
 259 }
 260
 261 //________________________________________________________________________
 262
 263 Double_t AliAnalysisTaskEmcalJetHMEC:: RelativePhi(Double_t mphi,Double_t vphi) {
 264
 265   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
 266   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
 267   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
 268   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
 269   Double_t dphi = vphi-mphi;
 270   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
 271   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
 272
 273   return dphi;//dphi in [-Pi, Pi]
 274 }
 275
 276 //________________________________________________________________________
 277 Int_t AliAnalysisTaskEmcalJetHMEC::GetCentBin(Double_t cent) const {
 278   // Get centrality bin.
 279
 280   Int_t centbin = -1;
 281   if (cent>=0 && cent<10) centbin = 0;
 282   else if (cent>=10 && cent<20) centbin = 1;
 283   else if (cent>=20 && cent<30) centbin = 2;
 284   else if (cent>=30 && cent<40) centbin = 3;
 285   else if (cent>=40 && cent<50) centbin = 4;
 286   else if (cent>=50 && cent<90) centbin = 5;
 287   return centbin;
 288 }
 289
 290 //________________________________________________________________________
 291 Int_t AliAnalysisTaskEmcalJetHMEC::GetEtaBin(Double_t eta) const {
 292   // Get eta bin for histos.
 293
 294   Int_t etabin = -1;
 295   if (TMath::Abs(eta)<=0.4) etabin = 0;
 296   else if (TMath::Abs(eta)>0.4 && TMath::Abs(eta)<0.8) etabin = 1;
 297   else if (TMath::Abs(eta)>=0.8) etabin = 2;
 298   return etabin;
 299 }
 300 //________________________________________________________________________
 301 Int_t AliAnalysisTaskEmcalJetHMEC::GetpTjetBin(Double_t pt) const
 302 {
 303   // Get jet pt  bin for histos.
 304
 305   Int_t ptbin = -1;
 306   if (pt>=15 && pt<20) ptbin = 0;
 307   else if (pt>=20 && pt<25) ptbin = 1;
 308   else if (pt>=25 && pt<30) ptbin = 2;
 309   else if (pt>=30 && pt<60) ptbin = 3;
 310   else if (pt>=60) ptbin = 4;
 311
 312   return ptbin;
 313 }
 314
 315 //________________________________________________________________________
 316 void AliAnalysisTaskEmcalJetHMEC::ExecOnce() {
 317   AliAnalysisTaskEmcalJet::ExecOnce();
 318
 319 }
 320
 321 //________________________________________________________________________
 322 Bool_t AliAnalysisTaskEmcalJetHMEC::Run() {
 323  // Main loop called for each event
 324   if(!fTracks){
 325     AliError(Form("No fTracks object!!\n"));
 326     return kTRUE;
 327   }
 328   if(!fJets){
 329     AliError(Form("No fJets object!!\n"));
 330     return kTRUE;
 331   }
 332
 333   // what kind of event do we have: AOD or ESD?
 334   Bool_t esdMode = kTRUE;
 335   if (dynamic_cast<AliAODEvent*>(InputEvent())) esdMode = kFALSE;
 336
 337   // if we have ESD event, set up ESD object
 338   if(esdMode){
 339     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 340     if (!fESD) {
 341       AliError(Form("ERROR: fESD not available\n"));
 342       return kTRUE;
 343     }
 344   }
 345
 346   // if we have AOD event, set up AOD object
 347   if(!esdMode){
 348     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 349     if(!fAOD) {
 350       AliError(Form("ERROR: fAOD not available\n"));
 351       return kTRUE;
 352     }
 353   }
 354
 355   TList *list = InputEvent()->GetList();
 356   if(!list) {
 357     AliError(Form("ERROR: list not attached\n"));
 358     return kTRUE;
 359   }
 360
 361   // get centrality
 362   if (fCent<0) {
 363     AliError(Form("Centrality negative: %f", fCent));
 364     return kTRUE;
 365   }
 366
 367   Int_t centbin = GetCentBin(fCent);
 368   if(centbin<0) return kTRUE;
 369
 370   Double_t fvertex[3]={0,0,0};
 371   InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
 372   Double_t zVtx=fvertex[2];
 373
 374   if(fabs(zVtx)>10.0) return kTRUE;
 375
 376   fHistCentrality->Fill(fCent);
 377
 378   TClonesArray *jets = 0;
 379   TClonesArray *tracks = 0;
 380
 381   tracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracks));
 382   if (!tracks) {
 383     AliError(Form("Pointer to tracks %s == 0", fTracks->GetName() ));
 384     return kTRUE;
 385   }
 386   const Int_t Ntracks=tracks->GetEntries();
 387
 388   jets= dynamic_cast<TClonesArray*>(list->FindObject(fJets));
 389   if (!jets) {
 390     AliError(Form("Pointer to tracks %s == 0", fJets->GetName() ));
 391     return kTRUE;
 392   }
 393   const Int_t Njets = jets->GetEntries();
 394
 395   //Leticia's loop to find hardest track
 396   Int_t iTT=-1;
 397   Double_t ptmax=-10;
 398
 399   for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++) {
 400     AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 401     if (!track) {
 402       printf("ERROR: Could not receive track %d\n", iTracks);
 403       continue;
 404     }
 405
 406     if(TMath::Abs(track->Eta())>0.9) continue;
 407     if(track->Pt()<0.15) continue;
 408     //iCount++;
 409     if(track->Pt()>ptmax){
 410           ptmax=track->Pt();
 411           iTT=iTracks;
 412     }
 413   }
 414
 415   Int_t ijethi=-1;
 416   Double_t highestjetpt=0.0;
 417   Int_t passedTTcut=0;
 418
 419   for (Int_t ijets = 0; ijets < Njets; ijets++){
 420       AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijets));
 421       if (!jet) continue;
 422       if(!AcceptthisJet(jet)) continue;
 423
 424       Double_t jetPt = jet->Pt();
 425
 426       if(highestjetpt<jetPt){
 427             ijethi=ijets;
 428             highestjetpt=jetPt;
 429       }
 430   }
 431
 432   // see if event is selected
 433   UInt_t trig = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 434
 435   for (Int_t ijet = 0; ijet < Njets; ijet++){
 436     AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
 437     if (!jet) continue;
 438
 439     // see if event is selected and our jet came from trigger event
 440     if (!(trig & fTriggerEventType)) continue;
 441     if (jet->Pt()<0.1) continue;
 442
 443     if(!AcceptthisJet(jet)) continue;
 444
 445     Double_t jetphi = jet->Phi();
 446     Double_t jetPt = jet->Pt();
 447     Double_t jeteta=jet->Eta();
 448
 449     Double_t leadjet=0;
 450     if (ijet==ijethi) leadjet=1;
 451
 452     fHistJetPt[centbin]->Fill(jet->Pt());
 453     fHistLeadJetPt[centbin]->Fill(jet->Pt());
 454
 455     if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
 456       fHistJetPtBias[centbin]->Fill(jet->Pt());
 457       fHistLeadJetPtBias[centbin]->Fill(jet->Pt());
 458     }
 459
 460     fHistJetEtaPhi->Fill(jet->Eta(),jetphi);
 461
 462     if(iTT>0){
 463           AliVTrack* TT = static_cast<AliVTrack*>(tracks->At(iTT));
 464           if(TMath::Abs(jetphi-TT->Phi()-TMath::Pi())<0.6) passedTTcut=1;
 465           else passedTTcut=0;
 466     }
 467
 468     if(passedTTcut)
 469           fHistJetPtTT[centbin]->Fill(jet->Pt());
 470
 471       Int_t iptjet=-1;
 472       iptjet=GetpTjetBin(jetPt);
 473       if(iptjet<0) continue;
 474
 475     if (jetPt > 15) {
 476
 477       for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++) {
 478             AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 479             if (!track) {
 480               printf("ERROR: Could not receive track %d\n", iTracks);
 481               continue;
 482             }
 483
 484             if(TMath::Abs(track->Eta())>fTrkEta) continue;
 485
 486             fHistTrackPt->Fill(track->Pt());
 487
 488             if (track->Pt()<0.15) continue;
 489
 490             Double_t trackphi = track->Phi();
 491             if (trackphi > TMath::Pi())
 492               trackphi = trackphi-2*TMath::Pi();
 493
 494             Double_t tracketa=track->Eta();
 495             Double_t trackpt=track->Pt();
 496             Double_t deta=tracketa-jeteta;
 497             Int_t ieta=GetEtaBin(deta);
 498             if (ieta<0) {
 499               AliError(Form("Eta Bin negative: %f", deta));
 500               continue;
 501             }
 502
 503             //Jet pt, track pt, dPhi,deta,fCent
 504             Double_t dphijh = RelativePhi(jetphi,trackphi);
 505             if (dphijh < -0.5*TMath::Pi())
 506               dphijh+= 2*TMath::Pi();
 507             if (dphijh > 1.5*TMath::Pi())
 508           dphijh-=2.*TMath::Pi();
 509
 510             fHistJetH[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
 511             fHistJetHEtaPhi->Fill(deta,dphijh);
 512
 513             Double_t dR=sqrt(deta*deta+dphijh*dphijh);
 514
 515             if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
 516               fHistJetHBias[centbin][iptjet][ieta]->Fill(dphijh,trackpt);
 517
 518           if(fDoLessSparseAxes) { // check if we want all dimensions
 519             Double_t triggerEntries[6] = {fCent,jetPt,trackpt,deta,dphijh,leadjet};
 520             fhnJH->Fill(triggerEntries);
 521           } else {
 522                 Double_t triggerEntries[8] = {fCent,jetPt,trackpt,deta,dphijh,leadjet,0.0,dR};
 523             fhnJH->Fill(triggerEntries);
 524           }
 525             }
 526
 527             if(passedTTcut)
 528               fHistJetHTT[centbin][iptjet][ieta]->Fill(dphijh,trackpt);
 529
 530           } //track loop
 531     }//jet pt cut
 532   }//jet loop
 533
 534   //Prepare to do event mixing
 535
 536   // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
 537   TObjArray* tracksClone = 0x0;
 538
 539   if(fDoEventMixing>0){
 540
 541     // event mixing
 542
 543     // 1. First get an event pool corresponding in mult (cent) and
 544     //    zvertex to the current event. Once initialized, the pool
 545     //    should contain nMix (reduced) events. This routine does not
 546     //    pre-scan the chain. The first several events of every chain
 547     //    will be skipped until the needed pools are filled to the
 548     //    specified depth. If the pool categories are not too rare, this
 549     //    should not be a problem. If they are rare, you could lose
 550     //    statistics.
 551
 552     // 2. Collect the whole pool's content of tracks into one TObjArray
 553     //    (bgTracks), which is effectively a single background super-event.
 554
 555     // 3. The reduced and bgTracks arrays must both be passed into
 556     //    FillCorrelations(). Also nMix should be passed in, so a weight
 557     //    of 1./nMix can be applied.
 558
 559     UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 560     // if event was not selected (triggered) for any reseason (should never happen) then return
 561     if (trigger==0)  return kTRUE;
 562
 563     AliEventPool* pool = fPoolMgr->GetEventPool(fCent, zVtx);
 564
 565     if (!pool){
 566       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCent, zVtx));
 567       return kTRUE;
 568     }
 569
 570 if(trigger & fTriggerEventType) {
 571     //check for a trigger jet
 572     if (pool->IsReady() || pool->NTracksInPool() > fNMIXtracks || pool->GetCurrentNEvents() >= fNMIXevents) {
 573
 574         for (Int_t ijet = 0; ijet < Njets; ijet++){
 575           Double_t leadjet=0;
 576           if (ijet==ijethi) leadjet=1;
 577
 578           AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
 579       if (!jet) continue;
 580
 581           if(!AcceptthisJet(jet)) continue;
 582
 583           Double_t jetPt = jet->Pt();
 584           Double_t jetphi = jet->Phi();
 585           Double_t jeteta=jet->Eta();
 586
 587       // make sure event contains jet above our threshold (reduce stats of sparse)
 588       if (jetPt < 15) continue;
 589
 590           Int_t nMix = pool->GetCurrentNEvents();
 591
 592           //Fill for biased jet triggers only
 593           if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
 594
 595             // Fill mixed-event histos here
 596             for (Int_t jMix=0; jMix<nMix; jMix++) {
 597                   TObjArray* bgTracks = pool->GetEvent(jMix);
 598                   const Int_t Nbgtrks = bgTracks->GetEntries();
 599
 600           for(Int_t ibg=0; ibg<Nbgtrks; ibg++){
 601                     AliPicoTrack *part = static_cast<AliPicoTrack*>(bgTracks->At(ibg));
 602                     if(!part) continue;
 603             if(TMath::Abs(part->Eta())>0.9) continue;
 604             if(part->Pt()<0.15) continue;
 605
 606                     Double_t DEta = part->Eta()-jeteta;
 607                     Double_t DPhi = RelativePhi(jetphi,part->Phi());
 608
 609                     Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
 610                     if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
 611                     if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
 612             if(fDoLessSparseAxes) {  // check if we want all the axis filled
 613               Double_t triggerEntries[6] = {fCent,jetPt,part->Pt(),DEta,DPhi,leadjet};
 614               fhnMixedEvents->Fill(triggerEntries,1./nMix);
 615             } else {
 616                       Double_t triggerEntries[8] = {fCent,jetPt,part->Pt(),DEta,DPhi,leadjet,0.0,DR};
 617               fhnMixedEvents->Fill(triggerEntries,1./nMix);
 618                     }
 619                     }
 620               }
 621             }
 622           }
 623     }
 624 }
 625
 626     if(trigger & fMixingEventType) {
 627       tracksClone = CloneAndReduceTrackList(tracks);
 628
 629       //update pool if jet in event or not
 630       pool->UpdatePool(tracksClone);
 631     }
 632
 633   } // end of event mixing
 634
 635   return kTRUE;
 636 }
 637
 638 //________________________________________________________________________
 639 void AliAnalysisTaskEmcalJetHMEC::Terminate(Option_t *)
 640 {
 641   //just terminate
 642 }
 643
 644 //________________________________________________________________________
 645 Int_t AliAnalysisTaskEmcalJetHMEC::AcceptthisJet(AliEmcalJet *jet)
 646 {
 647   //applies all jet cuts except pt
 648   float jetphi = jet->Phi();
 649   if (jetphi>TMath::Pi())
 650     jetphi = jetphi-2*TMath::Pi();
 651
 652   if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax))
 653     return 0;
 654   if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax))
 655     return 0;
 656   if (jet->Area()<fAreacut)
 657     return 0;
 658   //prevents 0 area jets from sneaking by when area cut == 0
 659   if (jet->Area()==0)
 660     return 0;
 661   //exclude jets with extremely high pt tracks which are likely misreconstructed
 662   if(jet->MaxTrackPt()>100)
 663     return 0;
 664
 665   //passed all above cuts
 666   return 1;
 667 }
 668
 669 //________________________________________________________________________
 670
 671 THnSparse* AliAnalysisTaskEmcalJetHMEC::NewTHnSparseF(const char* name, UInt_t entries){
 672    // generate new THnSparseF, axes are defined in GetDimParams()
 673
 674    Int_t count = 0;
 675    UInt_t tmp = entries;
 676    while(tmp!=0){
 677       count++;
 678       tmp = tmp &~ -tmp;  // clear lowest bit
 679    }
 680
 681    TString hnTitle(name);
 682    const Int_t dim = count;
 683    Int_t nbins[dim];
 684    Double_t xmin[dim];
 685    Double_t xmax[dim];
 686
 687    Int_t i=0;
 688    Int_t c=0;
 689    while(c<dim && i<32){
 690       if(entries&(1<<i)){
 691
 692          TString label("");
 693          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
 694          hnTitle += Form(";%s",label.Data());
 695          c++;
 696       }
 697
 698       i++;
 699    }
 700    hnTitle += ";";
 701
 702    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
 703 }
 704
 705 void AliAnalysisTaskEmcalJetHMEC::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
 706 {
 707    // stores label and binning of axis for THnSparse
 708
 709    const Double_t pi = TMath::Pi();
 710
 711    switch(iEntry){
 712
 713    case 0:
 714       label = "V0 centrality (%)";
 715          nbins = 10;
 716          xmin = 0.;
 717          xmax = 100.;
 718          break;
 719
 720    case 1:
 721       label = "corrected jet pt";
 722       nbins = 20;
 723       xmin = 0.;
 724       xmax = 200.;
 725       break;
 726
 727    case 2:
 728       if(fDoWiderTrackBin) {
 729         label = "track pT";
 730         nbins = 40;
 731         xmin = 0.;
 732         xmax = 10.;
 733       } else {
 734         label = "track pT";
 735         nbins = 100;
 736         xmin = 0.;
 737         xmax = 10;
 738       }
 739       break;
 740
 741    case 3:
 742       label = "deltaEta";
 743       nbins = 24;
 744       xmin = -1.2;
 745       xmax = 1.2;
 746       break;
 747
 748    case 4:
 749       label = "deltaPhi";
 750       nbins = 72;
 751       xmin = -0.5*pi;
 752       xmax = 1.5*pi;
 753       break;
 754
 755     case 5:
 756       label = "leading jet";
 757       nbins = 3;
 758       xmin = -0.5;
 759       xmax = 2.5;
 760       break;
 761
 762     case 6:
 763       label = "trigger track";
 764       nbins =10;
 765       xmin = 0;
 766       xmax = 50;
 767       break;
 768
 769     case 7:
 770       label = "deltaR";
 771       nbins = 10;
 772       xmin = 0.;
 773       xmax = 5.0;
 774       break;
 775
 776     case 8:
 777       label = "leading track";
 778       nbins = 13;
 779       xmin = 0;
 780       xmax = 50;
 781       break;
 782    }
 783 }
 784
 785 //_________________________________________________
 786 // From CF event mixing code PhiCorrelations
 787 TObjArray* AliAnalysisTaskEmcalJetHMEC::CloneAndReduceTrackList(TObjArray* tracks)
 788 {
 789   // clones a track list by using AliPicoTrack which uses much less memory (used for event mixing)
 790
 791   TObjArray* tracksClone = new TObjArray;
 792   tracksClone->SetOwner(kTRUE);
 793
 794   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
 795   {
 796     AliVParticle* particle = (AliVParticle*) tracks->At(i);
 797     if(TMath::Abs(particle->Eta())>fTrkEta) continue;
 798     if(particle->Pt()<0.15) continue;
 799
 800     Double_t trackpt=particle->Pt();
 801
 802     Int_t hadbin=-1;
 803     if(trackpt<0.5) hadbin=0;
 804     else if(trackpt<1) hadbin=1;
 805     else if(trackpt<2) hadbin=2;
 806     else if(trackpt<3) hadbin=3;
 807     else if(trackpt<5) hadbin=4;
 808     else if(trackpt<8) hadbin=5;
 809     else if(trackpt<20) hadbin=6;
 810
 811     if(hadbin>-1) fHistTrackEtaPhi[hadbin]->Fill(particle->Eta(),particle->Phi());
 812
 813     tracksClone->Add(new AliPicoTrack(particle->Pt(), particle->Eta(), particle->Phi(), particle->Charge(), 0, 0, 0, 0));
 814   }
 815
 816   return tracksClone;
 817 }