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