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