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