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