PWGJE/AliAnalysisTaskAj.cxx

   1
   2 // ******************************************
   3 // This task searches for events with large dijet imbalance
   4 // and then looks to the jet structure of the b-t-b jets.
   5 // *******************************************
   6
   7
   8 /**************************************************************************
   9  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  10  *                                                                        *
  11  * Author: The ALICE Off-line Project.                                    *
  12  * Contributors are mentioned in the code where appropriate.              *
  13  *                                                                        *
  14  * Permission to use, copy, modify and distribute this software and its   *
  15  * documentation strictly for non-commercial purposes is hereby granted   *
  16  * without fee, provided that the above copyright notice appears in all   *
  17  * copies and that both the copyright notice and this permission notice   *
  18  * appear in the supporting documentation. The authors make no claims     *
  19  * about the suitability of this software for any purpose. It is          *
  20  * provided "as is" without express or implied warranty.                  *
  21  **************************************************************************/
  22
  23
  24 #include "TChain.h"
  25 #include "TTree.h"
  26 #include "TMath.h"
  27 #include "TH1F.h"
  28 #include "TH2F.h"
  29 #include "TH3F.h"
  30 #include "THnSparse.h"
  31 #include "TMatrixD.h"
  32 #include "TMatrixDSym.h"
  33 #include "TMatrixDSymEigen.h"
  34
  35 #include "TCanvas.h"
  36
  37 #include "AliLog.h"
  38
  39 #include "AliAnalysisTask.h"
  40 #include "AliAnalysisManager.h"
  41
  42 #include "AliVEvent.h"
  43 #include "AliESDEvent.h"
  44 #include "AliESDInputHandler.h"
  45 #include "AliCentrality.h"
  46 #include "AliAnalysisHelperJetTasks.h"
  47 #include "AliInputEventHandler.h"
  48 #include "AliAODJetEventBackground.h"
  49 #include "AliAnalysisTaskFastEmbedding.h"
  50 #include "AliAODEvent.h"
  51 #include "AliAODHandler.h"
  52 #include "AliAODJet.h"
  53
  54 #include "AliAnalysisTaskAj.h"
  55
  56 ClassImp(AliAnalysisTaskAj)
  57
  58 AliAnalysisTaskAj::AliAnalysisTaskAj() :
  59 AliAnalysisTaskSE(),
  60 fESD(0x0),
  61 fAODIn(0x0),
  62 fAODOut(0x0),
  63 fAODExtension(0x0),
  64 fBackgroundBranch(""),
  65 fNonStdFile(""),
  66 fIsPbPb(kTRUE),
  67 fOfflineTrgMask(AliVEvent::kAny),
  68 fMinContribVtx(1),
  69 fVtxZMin(-10.),
  70 fVtxZMax(10.),
  71 fEvtClassMin(0),
  72 fEvtClassMax(4),
  73 fFilterMask(0),
  74 fRadioFrac(0.2),
  75 fMinDist(0.1),
  76 fCentMin(0.),
  77 fCentMax(100.),
  78 fNInputTracksMin(0),
  79 fNInputTracksMax(-1),
  80 fAngStructCloseTracks(0),
  81 fCheckMethods(0),
  82 fDoEventMixing(0),
  83 fJetEtaMin(-.5),
  84 fJetEtaMax(.5),
  85 fNevents(0x0),
  86 fTindex(0x0),
  87 fJetPtMin(20.),
  88 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
  89 fJetPtFractionMin(0.5),
  90 fNMatchJets(4),
  91 fMatchMaxDist(0.8),
  92 fKeepJets(kFALSE),
  93 fkNbranches(2),
  94 fkEvtClasses(12),
  95 fOutputList(0x0),
  96 fbEvent(kTRUE),
  97 fHistEvtSelection(0x0),
  98 fh2Pt1Pt2C10(0x0),
  99 fh2Pt1Pt2C40(0x0),
 100 fh3LocalCoordinates(0x0),
 101 fh2Sum2pt20(0x0),
 102 fh2Sum4pt20(0x0),
 103 fh2Sum6pt20(0x0),
 104 fh2Sum8pt20(0x0),
 105 fh2Sum12pt20(0x0),
 106 fh2Sum2lpt20(0x0),
 107 fh2Sum4lpt20(0x0),
 108 fh2Sum6lpt20(0x0),
 109 fh2Sum8lpt20(0x0),
 110 fh2Sum12lpt20(0x0),
 111 fh2Sum2pt40(0x0),
 112 fh2Sum4pt40(0x0),
 113 fh2Sum6pt40(0x0),
 114 fh2Sum8pt40(0x0),
 115 fh2Sum12pt40(0x0),
 116 fh2Sum2lpt40(0x0),
 117 fh2Sum4lpt40(0x0),
 118 fh2Sum6lpt40(0x0),
 119 fh2Sum8lpt40(0x0),
 120 fh2Sum12lpt40(0x0),
 121 fhnDeltaR(0x0)
 122  {
 123    // default Constructor
 124
 125
 126    fJetBranchName[0] = "";
 127    fJetBranchName[1] = "";
 128
 129    fListJets[0] = new TList;
 130    fListJets[1] = new TList;
 131 }
 132
 133 AliAnalysisTaskAj::AliAnalysisTaskAj(const char *name) :
 134 AliAnalysisTaskSE(name),
 135 fESD(0x0),
 136 fAODIn(0x0),
 137 fAODOut(0x0),
 138 fAODExtension(0x0),
 139 fBackgroundBranch(""),
 140 fNonStdFile(""),
 141 fIsPbPb(kTRUE),
 142 fOfflineTrgMask(AliVEvent::kAny),
 143 fMinContribVtx(1),
 144 fVtxZMin(-10.),
 145 fVtxZMax(10.),
 146 fEvtClassMin(0),
 147 fEvtClassMax(4),
 148 fFilterMask(0),
 149 fRadioFrac(0.2),
 150 fMinDist(0.1),
 151 fCentMin(0.),
 152 fCentMax(100.),
 153 fNInputTracksMin(0),
 154 fNInputTracksMax(-1),
 155 fAngStructCloseTracks(0),
 156 fCheckMethods(0),
 157 fDoEventMixing(0),
 158 fJetEtaMin(-.5),
 159 fJetEtaMax(.5),
 160 fNevents(0x0),
 161 fTindex(0x0),
 162 fJetPtMin(20.),
 163 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
 164 fJetPtFractionMin(0.5),
 165 fNMatchJets(4),
 166 fMatchMaxDist(0.8),
 167 fKeepJets(kFALSE),
 168 fkNbranches(2),
 169 fkEvtClasses(12),
 170 fOutputList(0x0),
 171 fbEvent(kTRUE),
 172 fHistEvtSelection(0x0),
 173 fh2Pt1Pt2C10(0x0),
 174 fh2Pt1Pt2C40(0x0),
 175 fh3LocalCoordinates(0x0),
 176 fh2Sum2pt20(0x0),
 177 fh2Sum4pt20(0x0),
 178 fh2Sum6pt20(0x0),
 179 fh2Sum8pt20(0x0),
 180 fh2Sum12pt20(0x0),
 181 fh2Sum2lpt20(0x0),
 182 fh2Sum4lpt20(0x0),
 183 fh2Sum6lpt20(0x0),
 184 fh2Sum8lpt20(0x0),
 185 fh2Sum12lpt20(0x0),
 186 fh2Sum2pt40(0x0),
 187 fh2Sum4pt40(0x0),
 188 fh2Sum6pt40(0x0),
 189 fh2Sum8pt40(0x0),
 190 fh2Sum12pt40(0x0),
 191 fh2Sum2lpt40(0x0),
 192 fh2Sum4lpt40(0x0),
 193 fh2Sum6lpt40(0x0),
 194 fh2Sum8lpt40(0x0),
 195 fh2Sum12lpt40(0x0),
 196 fhnDeltaR(0x0)
 197  {
 198    // Constructor
 199
 200
 201
 202
 203
 204
 205
 206    fJetBranchName[0] = "";
 207    fJetBranchName[1] = "";
 208
 209    fListJets[0] = new TList;
 210    fListJets[1] = new TList;
 211
 212    DefineOutput(1, TList::Class());
 213 }
 214
 215 AliAnalysisTaskAj::~AliAnalysisTaskAj()
 216 {
 217    delete fListJets[0];
 218    delete fListJets[1];
 219 }
 220
 221 void AliAnalysisTaskAj::SetBranchNames(const TString &branch1, const TString &branch2)
 222 {
 223    fJetBranchName[0] = branch1;
 224    fJetBranchName[1] = branch2;
 225 }
 226
 227 void AliAnalysisTaskAj::Init()
 228 {
 229
 230    // check for jet branches
 231    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 232       AliError("Jet branch name not set.");
 233    }
 234
 235 }
 236
 237 void AliAnalysisTaskAj::UserCreateOutputObjects()
 238 {
 239    // Create histograms
 240    // Called once
 241    OpenFile(1);
 242    if(!fOutputList) fOutputList = new TList;
 243    fOutputList->SetOwner(kTRUE);
 244
 245    Bool_t oldStatus = TH1::AddDirectoryStatus();
 246    TH1::AddDirectory(kFALSE);
 247
 248
 249    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
 250    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 251    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 252    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 253    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 254    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 255    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
 256
 257      UInt_t entries = 0; // bit coded, see GetDimParams() below
 258      entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7|1<<8;
 259      fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
 260
 261
 262      //change binning in centrality
 263      Double_t *xPt0 = new Double_t[4];
 264      xPt0[0] = 0.;
 265      for(int i = 1; i<=3;i++){
 266       if(xPt0[i-1]<10)xPt0[i] = xPt0[i-1] + 10; // 1 - 5
 267       else if(xPt0[i-1]<40)xPt0[i] = xPt0[i-1] + 30; // 5 - 12
 268       else xPt0[i] = xPt0[i-1] + 60.; // 18
 269      }
 270      fhnDeltaR->SetBinEdges(0,xPt0);
 271      delete [] xPt0;
 272
 273
 274
 275      //change binning in jTtrack
 276      Double_t *xPt3 = new Double_t[3];
 277      xPt3[0] = 0.;
 278      for(int i = 1; i<=2;i++){
 279       if(xPt3[i-1]<1)xPt3[i] = xPt3[i-1] + 1.; // 1 - 5
 280       else xPt3[i] = xPt3[i-1] + 149.; // 18
 281      }
 282      fhnDeltaR->SetBinEdges(3,xPt3);
 283      delete [] xPt3;
 284
 285
 286
 287      //change binning in pTtrack
 288      Double_t *xPt4 = new Double_t[5];
 289      xPt4[0] = 0.;
 290      for(int i = 1; i<=4;i++){
 291       if(xPt4[i-1]<0.4)xPt4[i] = xPt4[i-1] + 0.4; // 1 - 5
 292       else if(xPt4[i-1]<3)xPt4[i] = xPt4[i-1] + 2.6; // 5 - 12
 293       else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 7.4; // 5 - 12
 294       else xPt4[i] = xPt4[i-1] + 150.; // 18
 295      }
 296     fhnDeltaR->SetBinEdges(4,xPt4);
 297     delete [] xPt4;
 298
 299
 300
 301
 302
 303
 304     //change binning in HTI
 305      Double_t *xPt5 = new Double_t[4];
 306      xPt5[0] = 0.;
 307      for(int i = 1; i<=3;i++){
 308       if(xPt5[i-1]<10)xPt5[i] = xPt5[i-1] + 5; // 10 - 12
 309       else xPt5[i] = xPt5[i-1] + 40.; // 13
 310      }
 311     fhnDeltaR->SetBinEdges(8,xPt5);
 312     delete [] xPt5;
 313
 314
 315
 316     fh2Pt1Pt2C10 = new TH2F("Dijet spectra central","",20,0.,200.,20,0.,200.);
 317     fh2Pt1Pt2C40 = new TH2F("Dijet spectra peripheral","",20,0.,200.,20,0.,200.);
 318     fh3LocalCoordinates = new TH3F("Local coordinates","",10,-2,2,10,-2,2,10,0,100);
 319     fh2Sum2pt20 = new TH2F("pL R<0.2 pt20","",10,0.,1.,100,0.,100.);
 320     fh2Sum4pt20 = new TH2F("pL R<0.4 pt20","",10,0.,1.,100,0.,100.);
 321     fh2Sum6pt20 = new TH2F("pL R<0.6 pt20","",10,0.,1.,100,0.,100.);
 322     fh2Sum8pt20 = new TH2F("pL R<0.8 pt20","",10,0.,1.,100,0.,100.);
 323     fh2Sum12pt20 = new TH2F("pL R<1.2 pt20","",10,0.,1.,100,0.,100.);
 324     fh2Sum2lpt20 = new TH2F("pL R<0.2 low pt pt20","",10,0.,1.,100,0,100);
 325     fh2Sum4lpt20 = new TH2F("pL R<0.4 low pt pt20","",10,0.,1.,100,0,100);
 326     fh2Sum6lpt20 = new TH2F("pL R<0.6 low pt pt20","",10,0.,1.,100,0,100);
 327     fh2Sum8lpt20 = new TH2F("pL R<0.8 low pt pt20","",10,0.,1.,100,0,100);
 328     fh2Sum12lpt20 = new TH2F("pL R<1.2 low pt pt20","",10,0.,1.,100,0,100);
 329     fh2Sum2pt40 = new TH2F("pL R<0.2 pt40","",10,0.,1.,100,0.,100.);
 330     fh2Sum4pt40 = new TH2F("pL R<0.4 pt40","",10,0.,1.,100,0.,100.);
 331     fh2Sum6pt40 = new TH2F("pL R<0.6 pt40","",10,0.,1.,100,0.,100.);
 332     fh2Sum8pt40 = new TH2F("pL R<0.8 pt40","",10,0.,1.,100,0.,100.);
 333     fh2Sum12pt40 = new TH2F("pL R<1.2 pt40","",10,0.,1.,100,0.,100.);
 334     fh2Sum2lpt40 = new TH2F("pL R<0.2 low pt pt40","",10,0.,1.,100,0,100);
 335     fh2Sum4lpt40 = new TH2F("pL R<0.4 low pt pt40","",10,0.,1.,100,0,100);
 336     fh2Sum6lpt40 = new TH2F("pL R<0.6 low pt pt40","",10,0.,1.,100,0,100);
 337     fh2Sum8lpt40 = new TH2F("pL R<0.8 low pt pt40","",10,0.,1.,100,0,100);
 338     fh2Sum12lpt40 = new TH2F("pL R<1.2 low pt pt40","",10,0.,1.,100,0,100);
 339
 340
 341    fOutputList->Add(fHistEvtSelection);
 342    fOutputList->Add(fh2Pt1Pt2C10);
 343    fOutputList->Add(fh2Pt1Pt2C40);
 344    fOutputList->Add(fh3LocalCoordinates);
 345
 346    fOutputList->Add(fh2Sum2pt20);
 347    fOutputList->Add(fh2Sum4pt20);
 348    fOutputList->Add(fh2Sum6pt20);
 349    fOutputList->Add(fh2Sum8pt20);
 350    fOutputList->Add(fh2Sum12pt20);
 351    fOutputList->Add(fh2Sum2lpt20);
 352    fOutputList->Add(fh2Sum4lpt20);
 353    fOutputList->Add(fh2Sum6lpt20);
 354    fOutputList->Add(fh2Sum8lpt20);
 355    fOutputList->Add(fh2Sum12lpt20);
 356
 357    fOutputList->Add(fh2Sum2pt40);
 358    fOutputList->Add(fh2Sum4pt40);
 359    fOutputList->Add(fh2Sum6pt40);
 360    fOutputList->Add(fh2Sum8pt40);
 361    fOutputList->Add(fh2Sum12pt40);
 362    fOutputList->Add(fh2Sum2lpt40);
 363    fOutputList->Add(fh2Sum4lpt40);
 364    fOutputList->Add(fh2Sum6lpt40);
 365    fOutputList->Add(fh2Sum8lpt40);
 366    fOutputList->Add(fh2Sum12lpt40);
 367
 368    fOutputList->Add(fhnDeltaR);
 369
 370
 371
 372    // fOutputList->Add(fh3specbiased);
 373
 374    // =========== Switch on Sumw2 for all histos ===========
 375    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
 376       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 377       if (h1){
 378          h1->Sumw2();
 379          continue;
 380       }
 381     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
 382       if (hn){
 383          hn->Sumw2();
 384       }
 385    }
 386    TH1::AddDirectory(oldStatus);
 387
 388    PostData(1, fOutputList);
 389 }
 390
 391 void AliAnalysisTaskAj::UserExec(Option_t *)
 392 {
 393
 394
 395    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 396       AliError("Jet branch name not set.");
 397       return;
 398    }
 399
 400    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 401    if (!fESD) {
 402       AliError("ESD not available");
 403       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
 404       fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 405
 406       static AliAODEvent* aod = 0;
 407        // take all other information from the aod we take the tracks from
 408        if(!aod){
 409        if(!fESD)aod = fAODIn;
 410        else aod = fAODOut;}
 411
 412
 413
 414     if(fNonStdFile.Length()!=0){
 415     // case that we have an AOD extension we need can fetch the jets from the extended output
 416     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 417     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 418     if(!fAODExtension){
 419       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 420     }}
 421
 422
 423
 424
 425
 426    // -- event selection --
 427    fHistEvtSelection->Fill(1); // number of events before event selection
 428
 429    // physics selection
 430    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 431    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 432    cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
 433    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 434       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 435       fHistEvtSelection->Fill(2);
 436       PostData(1, fOutputList);
 437       return;
 438    }
 439
 440    // vertex selection
 441    if(!aod){
 442      if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
 443      fHistEvtSelection->Fill(3);
 444       PostData(1, fOutputList);
 445    }
 446    AliAODVertex* primVtx = aod->GetPrimaryVertex();
 447
 448    if(!primVtx){
 449      if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
 450      fHistEvtSelection->Fill(3);
 451      PostData(1, fOutputList);
 452      return;
 453    }
 454
 455    Int_t nTracksPrim = primVtx->GetNContributors();
 456    if ((nTracksPrim < fMinContribVtx) ||
 457          (primVtx->GetZ() < fVtxZMin) ||
 458          (primVtx->GetZ() > fVtxZMax) ){
 459       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 460       fHistEvtSelection->Fill(3);
 461       PostData(1, fOutputList);
 462       return;
 463    }
 464
 465    // event class selection (from jet helper task)
 466    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 467    if(fDebug) Printf("Event class %d", eventClass);
 468    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 469       fHistEvtSelection->Fill(4);
 470       PostData(1, fOutputList);
 471       return;
 472    }
 473
 474    // centrality selection
 475    AliCentrality *cent = 0x0;
 476    Double_t centValue = 0.;
 477    if(fESD) {cent = fESD->GetCentrality();
 478      if(cent) centValue = cent->GetCentralityPercentile("V0M");}
 479    else     centValue=aod->GetHeader()->GetCentrality();
 480
 481    if(fDebug) printf("centrality: %f\n", centValue);
 482    //  if (centValue < fCentMin || centValue > fCentMax){
 483    //  fHistEvtSelection->Fill(4);
 484    //  PostData(1, fOutputList);
 485    //  return;
 486    // }
 487
 488
 489    fHistEvtSelection->Fill(0);
 490    // accepted events
 491    // -- end event selection --
 492
 493    // get background
 494    AliAODJetEventBackground* externalBackground = 0;
 495    if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
 496       externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
 497       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 498    }
 499    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
 500      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
 501       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 502    }
 503
 504     if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
 505       externalBackground =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
 506       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 507     }
 508   Float_t rho = 0;
 509    if(externalBackground)rho = externalBackground->GetBackground(0);
 510
 511
 512    // fetch jets
 513    TClonesArray *aodJets[2];
 514    aodJets[0]=0;
 515    if(fAODOut&&!aodJets[0]){
 516    aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
 517    aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));  }
 518    if(fAODExtension && !aodJets[0]){
 519    aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
 520    aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
 521    if(fAODIn&&!aodJets[0]){
 522    aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
 523    aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
 524
 525    //Double_t ptsub[aodJets[0]->GetEntriesFast()];
 526    //Int_t inord[aodJets[0]->GetEntriesFast()];
 527    //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
 528    //  ptsub[n]=0;
 529    //  inord[n]=0;}
 530
 531    TList ParticleList;
 532    Int_t nT = GetListOfTracks(&ParticleList);
 533      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
 534       fListJets[iJetType]->Clear();
 535       if (!aodJets[iJetType]) continue;
 536
 537       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
 538
 539
 540       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 541          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 542          if (jet) fListJets[iJetType]->Add(jet);
 543          // if(iJetType==0){
 544          // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
 545       }}
 546
 547
 548    Double_t eta1=0;
 549    Int_t selec=-1;
 550    Double_t ptmax=-10;
 551    Double_t areaj=0;
 552    Double_t phij=0;
 553    Double_t etaj=0;
 554    Double_t ptj=0;
 555    Double_t ptcorrj=0;
 556    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
 557            AliAODJet* jetj = (AliAODJet*)(fListJets[0]->At(i));
 558            etaj  = jetj->Eta();
 559            phij  = jetj->Phi();
 560            ptj = jetj->Pt();
 561            if(ptj==0) continue;
 562            areaj = jetj->EffectiveAreaCharged();
 563            ptcorrj=ptj-rho*areaj;
 564            if(ptcorrj<=0) continue;
 565            if((etaj<fJetEtaMin)||(eta1>fJetEtaMax)) continue;
 566            if(ptcorrj>ptmax){ptmax=ptcorrj;
 567                              selec=i;}}
 568    ///hardest jet selected
 569     if(selec<0){PostData(1, fOutputList);
 570                return;}
 571     AliAODJet* jet1 = (AliAODJet*)(fListJets[0]->At(selec));
 572     //What is the hardest constituent track?
 573                        AliAODTrack* leadtrack1;
 574                        Int_t ippt=0;
 575                        Double_t ppt=-10;
 576                        TRefArray *genTrackList = jet1->GetRefTracks();
 577                        Int_t nTracksGenJet = genTrackList->GetEntriesFast();
 578                        AliAODTrack* genTrack;
 579                        for(Int_t ir=0; ir<nTracksGenJet; ++ir){
 580                        genTrack = (AliAODTrack*)(genTrackList->At(ir));
 581                        if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
 582                          ippt=ir;}}
 583                         leadtrack1=(AliAODTrack*)(genTrackList->At(ippt));
 584                         //If it doesn't exist or if it is greater that 100 GeV, discard.
 585                         if(!leadtrack1)  {PostData(1, fOutputList);
 586                           return;}
 587                         if(leadtrack1->Pt()>=100.){ PostData(1, fOutputList);
 588                           return;}
 589
 590            //Look to the back-to-back jet
 591            Int_t btb=-1;
 592            for(Int_t j=1;j<fListJets[0]->GetEntries();j++){
 593            if(j==selec) continue;
 594            AliAODJet* jetb = (AliAODJet*)(fListJets[0]->At(j));
 595            etaj  = jetb->Eta();
 596            phij  = jetb->Phi();
 597            ptj   = jetb->Pt();
 598            if(ptj<=0) continue;
 599            areaj = jetb->EffectiveAreaCharged();
 600            ptcorrj=ptj-rho*areaj;
 601            if(ptcorrj<=0) continue;
 602            if((etaj<fJetEtaMin)||(etaj>fJetEtaMax)) continue;
 603            Double_t dphij=RelativePhi(jetb->Phi(),jet1->Phi());
 604            if(TMath::Abs(dphij)>TMath::Pi()-0.2) { btb=j;
 605              break;}}
 606
 607             AliAODJet* jet2 = (AliAODJet*)(fListJets[0]->At(btb));
 608            //the back-to-back jet is also identified
 609
 610            if(btb<0){PostData(1, fOutputList);
 611                           return;}
 612
 613             Double_t ptcorr1=jet1->Pt()-rho*jet1->EffectiveAreaCharged();
 614             Double_t ptcorr2=jet2->Pt()-rho*jet2->EffectiveAreaCharged();
 615
 616             if(centValue<10.) fh2Pt1Pt2C10->Fill(ptcorr1,ptcorr2);
 617             if(centValue>40.) fh2Pt1Pt2C40->Fill(ptcorr1,ptcorr2);
 618
 619
 620            Double_t px2=jet2->Px();
 621            Double_t py2=jet2->Py();
 622            Double_t pz2=jet2->Pz();
 623            Double_t phi2=jet2->Phi();
 624            Double_t eta2=jet2->Eta();
 625            //Once we have have a dijet event,look to the structure of the back-to-back jet:
 626
 627       TVector3  ppJ1(px2, py2, pz2);
 628       TVector3  ppJ3(- px2 * pz2, - py2 * pz2, px2 * px2 + py2 * py2);
 629       ppJ3.SetMag(1.);
 630       TVector3  ppJ2(-py2, px2, 0);
 631       ppJ2.SetMag(1.);
 632       Float_t mxx    = 0.;
 633       Float_t myy    = 0.;
 634       Float_t mxy    = 0.;
 635       Int_t   nc     = 0;
 636       Float_t sump2  = 0.;
 637       Float_t ptMax  = 0.;
 638       Float_t etaMax = 0.;
 639       Float_t phiMax = 0.;
 640       Int_t   iMax   = -1;
 641       //1st loop over all tracks
 642          for(int it = 0;it<nT;++it){
 643          AliVParticle *track = (AliVParticle*)ParticleList.At(it);
 644          TVector3 pp(track->Px(), track->Py(), track->Pz());
 645           Float_t phi = track->Phi();
 646           Float_t eta = track->Eta();
 647           Float_t pt  = track->Pt();
 648           Float_t jT  = pp.Perp(ppJ1);
 649           Float_t deta = eta - eta2;
 650           Float_t dphi = phi - phi2;
 651           if (dphi >   TMath::Pi()) dphi =   2. * TMath::Pi() - dphi;
 652           if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
 653           Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
 654           /////////To compute the TM axis, we use particles with large jT
 655
 656               //cout<<"Jt spectrum............"<<ptbig<<" "<<jT<<endl;
 657           /////////We compute the TM with large jT particles only
 658              if((jT >1.)&&(r<1.2)){
 659             //longitudinal and perpendicular component of the track pT in the
 660             //local frame
 661               TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
 662               TVector3 pPerp = pp - pLong;
 663               //projection onto the two perpendicular vectors defined above
 664               Float_t ppjX = pPerp.Dot(ppJ2);
 665               Float_t ppjY = pPerp.Dot(ppJ3);
 666               Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
 667               //components of the sphericity matrix
 668               mxx += (ppjX * ppjX / ppjT);
 669               myy += (ppjY * ppjY / ppjT);
 670               mxy += (ppjX * ppjY / ppjT);
 671               nc++;
 672               sump2 += ppjT;
 673               // max pt
 674               if (pt > ptMax) {
 675                   ptMax  = pt;
 676                   iMax   = it;
 677                   etaMax = deta;
 678                   phiMax = dphi;
 679               }
 680                 } // R < 0.4
 681          } // 1st Track Loop
 682
 683   // At this point we have mxx, myy, mxy
 684       if (nc == 0) return;
 685 // Shericity Matrix
 686       const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
 687       TMatrixDSym m0(2,ele);
 688 // Find eigenvectors
 689       TMatrixDSymEigen m(m0);
 690       TVectorD eval(2);
 691       TMatrixD evecm = m.GetEigenVectors();
 692       eval  = m.GetEigenValues();
 693 // Largest eigenvector
 694       Int_t jev = 0;
 695       if (eval[0] < eval[1]) jev = 1;
 696       TVectorD evec0(2);
 697 // Principle axis
 698       evec0 = TMatrixDColumn(evecm, jev);
 699       TVector2 evec(evec0[0], evec0[1]);
 700 // Principle axis from leading partice
 701       Float_t phiM = TMath::ATan2(phiMax, etaMax);
 702       TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM));
 703       Float_t phistM = evecM.DeltaPhi(evec);
 704       if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
 705
 706       //////we have now the direction
 707       /////along which the sum of the projections of the particle
 708       ///momentum is higher.
 709       Double_t sum2lpt20=0;
 710       Double_t sum4lpt20=0;
 711       Double_t sum6lpt20=0;
 712       Double_t sum8lpt20=0;
 713       Double_t sum12lpt20=0;
 714       Double_t sum2pt20=0;
 715       Double_t sum4pt20=0;
 716       Double_t sum6pt20=0;
 717       Double_t sum8pt20=0;
 718       Double_t sum12pt20=0;
 719
 720        Double_t sum2lpt40=0;
 721       Double_t sum4lpt40=0;
 722       Double_t sum6lpt40=0;
 723       Double_t sum8lpt40=0;
 724       Double_t sum12lpt40=0;
 725       Double_t sum2pt40=0;
 726       Double_t sum4pt40=0;
 727       Double_t sum6pt40=0;
 728       Double_t sum8pt40=0;
 729       Double_t sum12pt40=0;
 730
 731           for (Int_t ip = 0; ip < nT; ip++) {
 732           AliVParticle *track = (AliVParticle*)ParticleList.At(ip);
 733           TVector3 pp(track->Px(), track->Py(), track->Pz());
 734           Float_t phi = track->Phi();
 735           Float_t eta = track->Eta();
 736           Float_t pt  = track->Pt();
 737           Float_t jT  = pp.Perp(ppJ1);
 738           Double_t deta = eta - eta2;
 739           Double_t deltaPhi = phi - phi2;
 740           Double_t r = TMath::Sqrt(deltaPhi * deltaPhi + deta * deta);
 741
 742             if(ptcorr2>20. && ptcorr2<40.){
 743             if(pt<0.4){
 744                  if(r<0.2) sum2lpt20=sum2lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());
 745                  if(r<0.4) sum4lpt20=sum4lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());
 746                  if(r<0.6) sum6lpt20=sum6lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());
 747                  if(r<0.8) sum8lpt20=sum8lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());
 748                  if(r<1.2) sum12lpt20=sum12lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());}
 749
 750                  if(r<0.2) sum2pt20=sum2pt20-1.*pt*TMath::Cos(phi-jet1->Phi());
 751                  if(r<0.4) sum4pt20=sum4pt20-1.*pt*TMath::Cos(phi-jet1->Phi());
 752                  if(r<0.6) sum6pt20=sum6pt20-1.*pt*TMath::Cos(phi-jet1->Phi());
 753                  if(r<0.8) sum8pt20=sum8pt20-1.*pt*TMath::Cos(phi-jet1->Phi());
 754                  if(r<1.2) sum12pt20=sum12pt20-1.*pt*TMath::Cos(phi-jet1->Phi());}
 755
 756
 757                 if(ptcorr2>40. && ptcorr2<60.){
 758             if(pt<0.4){
 759                  if(r<0.2) sum2lpt40=sum2lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());
 760                  if(r<0.4) sum4lpt40=sum4lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());
 761                  if(r<0.6) sum6lpt40=sum6lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());
 762                  if(r<0.8) sum8lpt40=sum8lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());
 763                  if(r<1.2) sum12lpt40=sum12lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());}
 764
 765                  if(r<0.2) sum2pt40=sum2pt40-1.*pt*TMath::Cos(phi-jet1->Phi());
 766                  if(r<0.4) sum4pt40=sum4pt40-1.*pt*TMath::Cos(phi-jet1->Phi());
 767                  if(r<0.6) sum6pt40=sum6pt40-1.*pt*TMath::Cos(phi-jet1->Phi());
 768                  if(r<0.8) sum8pt40=sum8pt40-1.*pt*TMath::Cos(phi-jet1->Phi());
 769                  if(r<1.2) sum12pt40=sum12pt40-1.*pt*TMath::Cos(phi-jet1->Phi());}
 770
 771
 772
 773
 774           TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
 775           TVector3 pPerp = pp - pLong;
 776           Float_t ppjX = pPerp.Dot(ppJ2);
 777           Float_t ppjY = pPerp.Dot(ppJ3);
 778           TVector2 vr(ppjX, ppjY) ;
 779           //and this is the angle between the particle and the TM axis.
 780           Float_t phistr = evec.DeltaPhi(vr);
 781
 782           //Double_t phistr=vr.Phi()-evec.Phi();
 783
 784           if(centValue<10.) fh3LocalCoordinates->Fill(ppjX,ppjY,ptcorr2);
 785           Double_t deltaEta = eta2-track->Eta();
 786
 787
 788           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
 789           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
 790           Double_t jetEntries[9] = {centValue,ptcorr1,ptcorr2,jT,pt,deltaEta,deltaPhi,phistr,ptMax};
 791           fhnDeltaR->Fill(jetEntries);
 792           }
 793           if(centValue<10.){
 794             if(ptcorr2>20.){
 795               if(ptcorr1>80.){
 796           Double_t aj=(ptcorr1-ptcorr2)/(ptcorr1+ptcorr2);
 797
 798           fh2Sum2pt20->Fill(aj,sum2pt20);
 799           fh2Sum4pt20->Fill(aj,sum4pt20);
 800           fh2Sum6pt20->Fill(aj,sum6pt20);
 801           fh2Sum8pt20->Fill(aj,sum8pt20);
 802           fh2Sum12pt20->Fill(aj,sum12pt20);
 803           fh2Sum2lpt20->Fill(aj,sum2lpt20);
 804           fh2Sum4lpt20->Fill(aj,sum4lpt20);
 805           fh2Sum6lpt20->Fill(aj,sum6lpt20);
 806           fh2Sum8lpt20->Fill(aj,sum8lpt20);
 807           fh2Sum12lpt20->Fill(aj,sum12lpt20);
 808
 809           fh2Sum2pt40->Fill(aj,sum2pt40);
 810           fh2Sum4pt40->Fill(aj,sum4pt40);
 811           fh2Sum6pt40->Fill(aj,sum6pt40);
 812           fh2Sum8pt40->Fill(aj,sum8pt40);
 813           fh2Sum12pt40->Fill(aj,sum12pt40);
 814           fh2Sum2lpt40->Fill(aj,sum2lpt40);
 815           fh2Sum4lpt40->Fill(aj,sum4lpt40);
 816           fh2Sum6lpt40->Fill(aj,sum6lpt40);
 817           fh2Sum8lpt40->Fill(aj,sum8lpt40);
 818           fh2Sum12lpt40->Fill(aj,sum12lpt40);
 819               }}}
 820
 821
 822
 823
 824
 825
 826
 827    PostData(1, fOutputList);
 828 }
 829
 830 void AliAnalysisTaskAj::Terminate(const Option_t *)
 831 {
 832    // Draw result to the screen
 833    // Called once at the end of the query
 834
 835    if (!GetOutputData(1))
 836    return;
 837 }
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849 Int_t  AliAnalysisTaskAj::GetListOfTracks(TList *list){
 850
 851     Int_t iCount = 0;
 852
 853      AliAODEvent *aod = 0;
 854      if(!fESD)aod = fAODIn;
 855      else aod = fAODOut;
 856      if(!aod) return iCount;
 857
 858     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 859       AliAODTrack *tr = aod->GetTrack(it);
 860       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 861       if(TMath::Abs(tr->Eta())>0.9)continue;
 862       if(tr->Pt()<0.15)continue;
 863       list->Add(tr);
 864       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
 865       iCount++;
 866     }
 867
 868
 869   return iCount;
 870
 871 }
 872
 873    Int_t  AliAnalysisTaskAj::GetHardestTrackBackToJet(AliAODJet *jetbig){
 874
 875       Int_t index=-1;
 876      AliAODEvent *aod = 0;
 877      if(!fESD)aod = fAODIn;
 878      else aod = fAODOut;
 879      if(!aod) return index;
 880
 881
 882     Double_t ptmax=-10;
 883     Double_t dphi=0;
 884     Double_t dif=0;
 885     Int_t iCount=0;
 886     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 887       AliAODTrack *tr = aod->GetTrack(it);
 888       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 889       if(TMath::Abs(tr->Eta())>0.9)continue;
 890       if(tr->Pt()<0.15)continue;
 891       iCount=iCount+1;
 892       dphi=RelativePhi(tr->Phi(),jetbig->Phi());
 893       if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
 894       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
 895       index=iCount-1;
 896       dif=dphi;  }}
 897
 898       return index;
 899
 900    }
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910  Int_t  AliAnalysisTaskAj::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
 911
 912     Int_t iCount = 0;
 913      AliAODEvent *aod = 0;
 914      if(!fESD)aod = fAODIn;
 915      else aod = fAODOut;
 916      if(!aod) return iCount;
 917
 918
 919
 920     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 921       AliAODTrack *tr = aod->GetTrack(it);
 922       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 923       if(TMath::Abs(tr->Eta())>0.9)continue;
 924       if(tr->Pt()<0.15)continue;
 925       Double_t disR=jetbig->DeltaR(tr);
 926       if(disR>0.8)  continue;
 927       list->Add(tr);
 928       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
 929       iCount++;
 930     }
 931
 932    list->Sort();
 933    return iCount;
 934
 935 }
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947 Int_t AliAnalysisTaskAj::GetNInputTracks()
 948 {
 949
 950    Int_t nInputTracks = 0;
 951       AliAODEvent *aod = 0;
 952      if(!fESD)aod = fAODIn;
 953      else aod = fAODOut;
 954      if(!aod) return nInputTracks;
 955
 956
 957    TString jbname(fJetBranchName[1]);
 958    //needs complete event, use jets without background subtraction
 959    for(Int_t i=1; i<=3; ++i){
 960       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
 961    }
 962    // use only HI event
 963    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
 964    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
 965
 966    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
 967    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
 968    if(!tmpAODjets){
 969       Printf("Jet branch %s not found", jbname.Data());
 970       Printf("AliAnalysisTaskAj::GetNInputTracks FAILED");
 971       return -1;
 972    }
 973
 974    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
 975       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
 976       if(!jet) continue;
 977       TRefArray *trackList = jet->GetRefTracks();
 978       Int_t nTracks = trackList->GetEntriesFast();
 979       nInputTracks += nTracks;
 980       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
 981    }
 982    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
 983
 984    return nInputTracks;
 985 }
 986
 987
 988
 989 Double_t AliAnalysisTaskAj::RelativePhi(Double_t mphi,Double_t vphi){
 990
 991   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
 992   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
 993   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
 994   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
 995   double dphi = mphi-vphi;
 996   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
 997   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
 998   return dphi;//dphi in [-Pi, Pi]
 999 }
1000
1001
1002
1003 THnSparse* AliAnalysisTaskAj::NewTHnSparseF(const char* name, UInt_t entries)
1004 {
1005    // generate new THnSparseF, axes are defined in GetDimParams()
1006
1007    Int_t count = 0;
1008    UInt_t tmp = entries;
1009    while(tmp!=0){
1010       count++;
1011       tmp = tmp &~ -tmp;  // clear lowest bit
1012    }
1013
1014    TString hnTitle(name);
1015    const Int_t dim = count;
1016    Int_t nbins[dim];
1017    Double_t xmin[dim];
1018    Double_t xmax[dim];
1019
1020    Int_t i=0;
1021    Int_t c=0;
1022    while(c<dim && i<32){
1023       if(entries&(1<<i)){
1024
1025          TString label("");
1026          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1027          hnTitle += Form(";%s",label.Data());
1028          c++;
1029       }
1030
1031       i++;
1032    }
1033    hnTitle += ";";
1034
1035    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1036 }
1037
1038 void AliAnalysisTaskAj::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1039 {
1040    // stores label and binning of axis for THnSparse
1041
1042    const Double_t pi = TMath::Pi();
1043
1044    switch(iEntry){
1045
1046    case 0:
1047       label = "V0 centrality (%)";
1048
1049          nbins = 3;
1050          xmin = 0.;
1051          xmax = 100.;
1052          break;
1053
1054
1055    case 1:
1056       label = "corrected jet pt1";
1057          nbins = 10;
1058          xmin = 0.;
1059          xmax = 200.;
1060           break;
1061
1062       case 2:
1063       label = "corrected jet pt2";
1064          nbins = 10;
1065          xmin = 0.;
1066          xmax = 200.;
1067           break;
1068
1069
1070
1071    case 3:
1072       label = "track jT";
1073
1074          nbins = 2;
1075          xmin = 0.;
1076          xmax = 150.;
1077          break;
1078
1079    case 4:
1080       label = "track pT";
1081
1082          nbins = 4;
1083          xmin = 0.;
1084          xmax = 150.;
1085          break;
1086
1087
1088       case 5:
1089       label = "deltaEta";
1090       nbins = 8;
1091       xmin = -1.6;
1092       xmax = 1.6;
1093       break;
1094
1095
1096       case 6:
1097       label = "deltaPhi";
1098       nbins = 60;
1099       xmin = -0.5*pi;
1100       xmax = 1.5*pi;
1101       break;
1102
1103       case 7:
1104       label = "deltaPhiTM";
1105       nbins = 60;
1106       xmin = 0.;
1107       xmax = 1.3*pi;
1108       break;
1109
1110
1111
1112       case 8:
1113       label = "leading track";
1114       nbins = 3;
1115       xmin = 0;
1116       xmax = 50;
1117       break;
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128    }
1129
1130 }
1131