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