PWGJE/AliAnalysisTaskJetCore.cxx

   1
   2 // ******************************************
   3 // This task computes several jet observables like
   4 // the fraction of energy in inner and outer coronnas,
   5 // jet-track correlations,triggered jet shapes and
   6 // correlation strength distribution of particles inside jets.
   7 // Author: lcunquei@cern.ch
   8 // *******************************************
   9
  10
  11 /**************************************************************************
  12  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  13  *                                                                        *
  14  * Author: The ALICE Off-line Project.                                    *
  15  * Contributors are mentioned in the code where appropriate.              *
  16  *                                                                        *
  17  * Permission to use, copy, modify and distribute this software and its   *
  18  * documentation strictly for non-commercial purposes is hereby granted   *
  19  * without fee, provided that the above copyright notice appears in all   *
  20  * copies and that both the copyright notice and this permission notice   *
  21  * appear in the supporting documentation. The authors make no claims     *
  22  * about the suitability of this software for any purpose. It is          *
  23  * provided "as is" without express or implied warranty.                  *
  24  **************************************************************************/
  25
  26
  27 #include "TChain.h"
  28 #include "TTree.h"
  29 #include "TMath.h"
  30 #include "TH1F.h"
  31 #include "TH2F.h"
  32 #include "TH3F.h"
  33 #include "THnSparse.h"
  34 #include "TCanvas.h"
  35
  36 #include "AliLog.h"
  37
  38 #include "AliAnalysisTask.h"
  39 #include "AliAnalysisManager.h"
  40
  41 #include "AliVEvent.h"
  42 #include "AliESDEvent.h"
  43 #include "AliESDInputHandler.h"
  44 #include "AliCentrality.h"
  45 #include "AliAnalysisHelperJetTasks.h"
  46 #include "AliInputEventHandler.h"
  47 #include "AliAODJetEventBackground.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliAnalysisTaskFastEmbedding.h"
  50 #include "AliAODEvent.h"
  51 #include "AliAODHandler.h"
  52 #include "AliAODJet.h"
  53
  54 #include "AliAnalysisTaskJetCore.h"
  55
  56 using std::cout;
  57 using std::endl;
  58
  59 ClassImp(AliAnalysisTaskJetCore)
  60
  61 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
  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 fFilterMaskBestPt(0),
  78 fFilterType(0),
  79 fRadioFrac(0.2),
  80 fMinDist(0.1),
  81 fCentMin(0.),
  82 fCentMax(100.),
  83 fNInputTracksMin(0),
  84 fNInputTracksMax(-1),
  85 fAngStructCloseTracks(0),
  86 fCheckMethods(0),
  87 fDoEventMixing(0),
  88 fFlagPhiBkg(0),
  89 fFlagEtaBkg(0),
  90 fFlagJetHadron(0),
  91 fFlagRandom(0),
  92 fFlagOnlyRecoil(0),
  93 fFlagOnlyHardest(1),
  94 fTrackTypeRec(kTrackUndef),
  95 fRPAngle(0),
  96 fNRPBins(50),
  97 fJetEtaMin(-.5),
  98 fJetEtaMax(.5),
  99 fNevents(0),
 100 fTindex(0),
 101 fTrigBufferIndex(0),
 102 fCountAgain(0),
 103 fJetPtMin(20.),
 104 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
 105 fJetPtFractionMin(0.5),
 106 fNMatchJets(4),
 107 fMatchMaxDist(0.8),
 108 fKeepJets(kFALSE),
 109 fkNbranches(2),
 110 fkEvtClasses(12),
 111 fOutputList(0x0),
 112 fbEvent(kTRUE),
 113 fHistEvtSelection(0x0),
 114 fhnDeltaR(0x0),
 115 fhnMixedEvents(0x0),
 116 fh2JetCoreMethod1C10(0x0),
 117 fh2JetCoreMethod2C10(0x0),
 118 fh2JetCoreMethod1C20(0x0),
 119 fh2JetCoreMethod2C20(0x0),
 120 fh2JetCoreMethod1C30(0x0),
 121 fh2JetCoreMethod2C30(0x0),
 122 fh2JetCoreMethod1C60(0x0),
 123 fh2JetCoreMethod2C60(0x0),
 124 fh3JetTrackC3060(0x0),
 125 fh3JetTrackC20(0x0),
 126 fh2AngStructpt1C10(0x0),
 127 fh2AngStructpt2C10(0x0),
 128 fh2AngStructpt3C10(0x0),
 129 fh2AngStructpt4C10(0x0),
 130 fh2AngStructpt1C20(0x0),
 131 fh2AngStructpt2C20(0x0),
 132 fh2AngStructpt3C20(0x0),
 133 fh2AngStructpt4C20(0x0),
 134 fh2AngStructpt1C30(0x0),
 135 fh2AngStructpt2C30(0x0),
 136 fh2AngStructpt3C30(0x0),
 137 fh2AngStructpt4C30(0x0),
 138 fh2AngStructpt1C60(0x0),
 139 fh2AngStructpt2C60(0x0),
 140 fh2AngStructpt3C60(0x0),
 141 fh2AngStructpt4C60(0x0),
 142 fh2Ntriggers(0x0),
 143 fh2Ntriggers2C10(0x0),
 144 fh2Ntriggers2C20(0x0),
 145 fh3JetDensity(0x0),
 146 fh3JetDensityA4(0x0),
 147 fh2RPJetsC10(0x0),
 148 fh2RPJetsC20(0x0),
 149 fh2RPTC10(0x0),
 150 fh2RPTC20(0x0),
 151 fh3spectriggeredC10(0x0),
 152 fh3spectriggeredC20(0x0),
 153 fh3spectriggeredC3060(0x0)
 154
 155
 156 {
 157    // default Constructor
 158
 159
 160  // Trigger buffer.
 161    for(Int_t i=0; i<10; i++) {
 162                 for(Int_t j=0; j<6; j++) {
 163                         fTrigBuffer[i][j]=0;
 164                 }
 165         }
 166
 167
 168
 169
 170
 171    fJetBranchName[0] = "";
 172    fJetBranchName[1] = "";
 173
 174    fListJets[0] = new TList;
 175    fListJets[1] = new TList;
 176 }
 177
 178 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
 179 AliAnalysisTaskSE(name),
 180 fESD(0x0),
 181 fAODIn(0x0),
 182 fAODOut(0x0),
 183 fAODExtension(0x0),
 184 fBackgroundBranch(""),
 185 fNonStdFile(""),
 186 fIsPbPb(kTRUE),
 187 fOfflineTrgMask(AliVEvent::kAny),
 188 fMinContribVtx(1),
 189 fVtxZMin(-10.),
 190 fVtxZMax(10.),
 191 fEvtClassMin(0),
 192 fEvtClassMax(4),
 193 fFilterMask(0),
 194 fFilterMaskBestPt(0),
 195 fFilterType(0),
 196 fRadioFrac(0.2),
 197 fMinDist(0.1),
 198 fCentMin(0.),
 199 fCentMax(100.),
 200 fNInputTracksMin(0),
 201 fNInputTracksMax(-1),
 202 fAngStructCloseTracks(0),
 203 fCheckMethods(0),
 204 fDoEventMixing(0),
 205 fFlagPhiBkg(0),
 206 fFlagEtaBkg(0),
 207 fFlagJetHadron(0),
 208 fFlagRandom(0),
 209 fFlagOnlyRecoil(0),
 210 fFlagOnlyHardest(1),
 211 fTrackTypeRec(kTrackUndef),
 212 fRPAngle(0),
 213 fNRPBins(50),
 214 fJetEtaMin(-.5),
 215 fJetEtaMax(.5),
 216 fNevents(0),
 217 fTindex(0),
 218 fTrigBufferIndex(0),
 219 fCountAgain(0),
 220 fJetPtMin(20.),
 221 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
 222 fJetPtFractionMin(0.5),
 223 fNMatchJets(4),
 224 fMatchMaxDist(0.8),
 225 fKeepJets(kFALSE),
 226 fkNbranches(2),
 227 fkEvtClasses(12),
 228 fOutputList(0x0),
 229 fbEvent(kTRUE),
 230 fHistEvtSelection(0x0),
 231 fhnDeltaR(0x0),
 232 fhnMixedEvents(0x0),
 233 fh2JetCoreMethod1C10(0x0),
 234 fh2JetCoreMethod2C10(0x0),
 235 fh2JetCoreMethod1C20(0x0),
 236 fh2JetCoreMethod2C20(0x0),
 237 fh2JetCoreMethod1C30(0x0),
 238 fh2JetCoreMethod2C30(0x0),
 239 fh2JetCoreMethod1C60(0x0),
 240 fh2JetCoreMethod2C60(0x0),
 241 fh3JetTrackC3060(0x0),
 242 fh3JetTrackC20(0x0),
 243 fh2AngStructpt1C10(0x0),
 244 fh2AngStructpt2C10(0x0),
 245 fh2AngStructpt3C10(0x0),
 246 fh2AngStructpt4C10(0x0),
 247 fh2AngStructpt1C20(0x0),
 248 fh2AngStructpt2C20(0x0),
 249 fh2AngStructpt3C20(0x0),
 250 fh2AngStructpt4C20(0x0),
 251 fh2AngStructpt1C30(0x0),
 252 fh2AngStructpt2C30(0x0),
 253 fh2AngStructpt3C30(0x0),
 254 fh2AngStructpt4C30(0x0),
 255 fh2AngStructpt1C60(0x0),
 256 fh2AngStructpt2C60(0x0),
 257 fh2AngStructpt3C60(0x0),
 258 fh2AngStructpt4C60(0x0),
 259 fh2Ntriggers(0x0),
 260 fh2Ntriggers2C10(0x0),
 261 fh2Ntriggers2C20(0x0),
 262 fh3JetDensity(0x0),
 263 fh3JetDensityA4(0x0),
 264 fh2RPJetsC10(0x0),
 265 fh2RPJetsC20(0x0),
 266 fh2RPTC10(0x0),
 267 fh2RPTC20(0x0),
 268 fh3spectriggeredC10(0x0),
 269 fh3spectriggeredC20(0x0),
 270 fh3spectriggeredC3060(0x0)
 271
 272  {
 273    // Constructor
 274
 275
 276     for(Int_t i=0; i<10; i++) {
 277        for(Int_t j=0; j<6; j++) {
 278             fTrigBuffer[i][j]=0;
 279                 }
 280     }
 281
 282
 283
 284    fJetBranchName[0] = "";
 285    fJetBranchName[1] = "";
 286
 287    fListJets[0] = new TList;
 288    fListJets[1] = new TList;
 289
 290    DefineOutput(1, TList::Class());
 291 }
 292
 293 AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
 294 {
 295    delete fListJets[0];
 296    delete fListJets[1];
 297 }
 298
 299 void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
 300 {
 301    fJetBranchName[0] = branch1;
 302    fJetBranchName[1] = branch2;
 303 }
 304
 305 void AliAnalysisTaskJetCore::Init()
 306 {
 307
 308    // check for jet branches
 309    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 310       AliError("Jet branch name not set.");
 311    }
 312
 313 }
 314
 315 void AliAnalysisTaskJetCore::UserCreateOutputObjects()
 316 {
 317    // Create histograms
 318    // Called once
 319    OpenFile(1);
 320    if(!fOutputList) fOutputList = new TList;
 321    fOutputList->SetOwner(kTRUE);
 322
 323    Bool_t oldStatus = TH1::AddDirectoryStatus();
 324    TH1::AddDirectory(kFALSE);
 325
 326
 327    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
 328    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 329    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 330    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 331    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 332    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 333    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
 334
 335      UInt_t entries = 0; // bit coded, see GetDimParams() below
 336      entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7;
 337      fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
 338
 339      //change binning in pTtrack
 340      Double_t *xPt3 = new Double_t[10];
 341      xPt3[0] = 0.;
 342      for(int i = 1; i<=9;i++){
 343       if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.4; // 1 - 5
 344       else if(xPt3[i-1]<11)xPt3[i] = xPt3[i-1] + 3; // 5 - 12
 345       else xPt3[i] = xPt3[i-1] + 150.; // 18
 346      }
 347     fhnDeltaR->SetBinEdges(2,xPt3);
 348     delete [] xPt3;
 349
 350     //change binning in HTI
 351      Double_t *xPt4 = new Double_t[14];
 352      xPt4[0] = 0.;
 353      for(int i = 1; i<=13;i++){
 354       if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 1; // 1 - 10
 355       else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] + 5; // 10 - 12
 356       else xPt4[i] = xPt4[i-1] + 30.; // 13
 357      }
 358     fhnDeltaR->SetBinEdges(6,xPt4);
 359     delete [] xPt4;
 360
 361
 362
 363
 364
 365
 366
 367      if(fDoEventMixing){
 368      UInt_t cifras = 0; // bit coded, see GetDimParams() below
 369      cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7;
 370      fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);}
 371
 372     if(fCheckMethods){
 373     fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
 374     fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
 375     fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
 376     fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
 377     fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
 378     fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
 379     fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
 380     fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);}
 381
 382     fh3JetTrackC3060=new TH3F("JetTrackC3060","",50,0,50,150,0.,150.,35,0.,3.5);
 383     fh3JetTrackC20=new TH3F("JetTrackC20","",50,0,50,150,0.,150.,35,0.,3.5);
 384     if(fAngStructCloseTracks>0){
 385     fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,10.);
 386     fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,10.);
 387     fh2AngStructpt3C10 = new TH2F("Ang struct pt3 C10","",15,0.,1.5,150,0.,10.);
 388     fh2AngStructpt4C10 = new TH2F("Ang struct pt4 C10","",15,0.,1.5,150,0.,10.);
 389     fh2AngStructpt1C20 = new TH2F("Ang struct pt1 C20","",15,0.,1.5,150,0.,10.);
 390     fh2AngStructpt2C20 = new TH2F("Ang struct pt2 C20","",15,0.,1.5,150,0.,10.);
 391     fh2AngStructpt3C20 = new TH2F("Ang struct pt3 C20","",15,0.,1.5,150,0.,10.);
 392     fh2AngStructpt4C20 = new TH2F("Ang struct pt4 C20","",15,0.,1.5,150,0.,10.);
 393     fh2AngStructpt1C30 = new TH2F("Ang struct pt1 C30","",15,0.,1.5,150,0.,10.);
 394     fh2AngStructpt2C30 = new TH2F("Ang struct pt2 C30","",15,0.,1.5,150,0.,10.);
 395     fh2AngStructpt3C30 = new TH2F("Ang struct pt3 C30","",15,0.,1.5,150,0.,10.);
 396     fh2AngStructpt4C30 = new TH2F("Ang struct pt4 C30","",15,0.,1.5,150,0.,10.);
 397     fh2AngStructpt1C60 = new TH2F("Ang struct pt1 C60","",15,0.,1.5,150,0.,10.);
 398     fh2AngStructpt2C60 = new TH2F("Ang struct pt2 C60","",15,0.,1.5,150,0.,10.);
 399     fh2AngStructpt3C60 = new TH2F("Ang struct pt3 C60","",15,0.,1.5,150,0.,10.);
 400     fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,10.); }
 401
 402
 403
 404     fh2Ntriggers=new TH2F("# of triggers","",10,0.,100.,50,0.,50.);
 405     fh2Ntriggers2C10=new TH2F("# of triggers2C10","",50,0.,50.,50,0.,50.);
 406     fh2Ntriggers2C20=new TH2F("# of triggers2C20","",50,0.,50.,50,0.,50.);
 407     fh3JetDensity=new TH3F("Jet density vs mutliplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
 408     fh3JetDensityA4=new TH3F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
 409     fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,35,100,0.,100.);
 410     fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,35,100,0.,100.);
 411     fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,35,50,0.,50.);
 412     fh2RPTC20=new TH2F("RPTriggerC20","",35,0.,35,50,0.,50.);
 413     fh3spectriggeredC10 = new TH3F("Triggered spectrumC10","",100,0.,1.,140,-80.,200.,50,0.,50.);
 414     fh3spectriggeredC20 = new TH3F("Triggered spectrumC20","",100,0.,1.,140,-80.,200.,50,0.,50.);
 415     fh3spectriggeredC3060 = new TH3F("Triggered spectrumC3060","",100,0.,1.,140,-80.,200.,10,0.,50.);
 416
 417
 418
 419    fOutputList->Add(fHistEvtSelection);
 420
 421    fOutputList->Add(fhnDeltaR);
 422
 423    fOutputList->Add(fhnMixedEvents);
 424
 425
 426
 427       if(fCheckMethods){
 428
 429       fOutputList->Add(fh2JetCoreMethod1C10);
 430       fOutputList->Add(fh2JetCoreMethod2C10);
 431       fOutputList->Add(fh2JetCoreMethod1C20);
 432       fOutputList->Add(fh2JetCoreMethod2C20);
 433       fOutputList->Add(fh2JetCoreMethod1C30);
 434       fOutputList->Add(fh2JetCoreMethod2C30);
 435       fOutputList->Add(fh2JetCoreMethod1C60);
 436       fOutputList->Add(fh2JetCoreMethod2C60);}
 437
 438       fOutputList->Add(fh3JetTrackC3060);
 439       fOutputList->Add(fh3JetTrackC20);
 440
 441
 442
 443
 444
 445         if(fAngStructCloseTracks>0){
 446        fOutputList->Add(fh2AngStructpt1C10);
 447        fOutputList->Add(fh2AngStructpt2C10);
 448        fOutputList->Add(fh2AngStructpt3C10);
 449        fOutputList->Add(fh2AngStructpt4C10);
 450        fOutputList->Add(fh2AngStructpt1C20);
 451        fOutputList->Add(fh2AngStructpt2C20);
 452        fOutputList->Add(fh2AngStructpt3C20);
 453        fOutputList->Add(fh2AngStructpt4C20);
 454        fOutputList->Add(fh2AngStructpt1C30);
 455        fOutputList->Add(fh2AngStructpt2C30);
 456        fOutputList->Add(fh2AngStructpt3C30);
 457        fOutputList->Add(fh2AngStructpt4C30);
 458        fOutputList->Add(fh2AngStructpt1C60);
 459        fOutputList->Add(fh2AngStructpt2C60);
 460        fOutputList->Add(fh2AngStructpt3C60);
 461        fOutputList->Add(fh2AngStructpt4C60);}
 462
 463
 464
 465
 466
 467         fOutputList->Add(fh2Ntriggers);
 468         fOutputList->Add(fh2Ntriggers2C10);
 469         fOutputList->Add(fh2Ntriggers2C20);
 470         fOutputList->Add(fh3JetDensity);
 471         fOutputList->Add(fh3JetDensityA4);
 472         fOutputList->Add(fh2RPJetsC10);
 473         fOutputList->Add(fh2RPJetsC20);
 474          fOutputList->Add(fh2RPTC10);
 475         fOutputList->Add(fh2RPTC20);
 476         fOutputList->Add(fh3spectriggeredC10);
 477         fOutputList->Add(fh3spectriggeredC20);
 478         fOutputList->Add(fh3spectriggeredC3060);
 479
 480
 481    // =========== Switch on Sumw2 for all histos ===========
 482    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
 483       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 484       if (h1){
 485          h1->Sumw2();
 486          continue;
 487       }
 488     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
 489       if (hn){
 490          hn->Sumw2();
 491       }
 492    }
 493    TH1::AddDirectory(oldStatus);
 494
 495    PostData(1, fOutputList);
 496 }
 497
 498 void AliAnalysisTaskJetCore::UserExec(Option_t *)
 499 {
 500
 501
 502    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 503       AliError("Jet branch name not set.");
 504       return;
 505    }
 506
 507    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 508    if (!fESD) {
 509       AliError("ESD not available");
 510       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 511    }
 512       fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 513
 514        static AliAODEvent* aod = 0;
 515        // take all other information from the aod we take the tracks from
 516        if(!aod){
 517        if(!fESD)aod = fAODIn;
 518        else aod = fAODOut;}
 519
 520
 521
 522     if(fNonStdFile.Length()!=0){
 523     // case that we have an AOD extension we need can fetch the jets from the extended output
 524     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 525     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 526     if(!fAODExtension){
 527       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 528     }}
 529
 530
 531    // -- event selection --
 532    fHistEvtSelection->Fill(1); // number of events before event selection
 533
 534    // physics selection
 535    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 536    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 537          std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
 538    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 539       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 540       fHistEvtSelection->Fill(2);
 541       PostData(1, fOutputList);
 542       return;
 543    }
 544
 545    // vertex selection
 546    if(!aod){
 547      if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
 548      fHistEvtSelection->Fill(3);
 549       PostData(1, fOutputList);
 550    }
 551    AliAODVertex* primVtx = aod->GetPrimaryVertex();
 552
 553    if(!primVtx){
 554      if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
 555      fHistEvtSelection->Fill(3);
 556      PostData(1, fOutputList);
 557      return;
 558    }
 559
 560    Int_t nTracksPrim = primVtx->GetNContributors();
 561    if ((nTracksPrim < fMinContribVtx) ||
 562          (primVtx->GetZ() < fVtxZMin) ||
 563          (primVtx->GetZ() > fVtxZMax) ){
 564       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 565       fHistEvtSelection->Fill(3);
 566       PostData(1, fOutputList);
 567       return;
 568    }
 569
 570    // event class selection (from jet helper task)
 571    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 572    if(fDebug) Printf("Event class %d", eventClass);
 573    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 574       fHistEvtSelection->Fill(4);
 575       PostData(1, fOutputList);
 576       return;
 577    }
 578
 579    // centrality selection
 580    AliCentrality *cent = 0x0;
 581    Double_t centValue = 0.;
 582    if(fIsPbPb){
 583    if(fESD) {cent = fESD->GetCentrality();
 584      if(cent) centValue = cent->GetCentralityPercentile("V0M");}
 585    else     centValue=aod->GetHeader()->GetCentrality();
 586
 587    if(fDebug) printf("centrality: %f\n", centValue);
 588       if (centValue < fCentMin || centValue > fCentMax){
 589       fHistEvtSelection->Fill(4);
 590       PostData(1, fOutputList);
 591       return;
 592       }}
 593
 594
 595    fHistEvtSelection->Fill(0);
 596    // accepted events
 597    // -- end event selection --
 598
 599    // get background
 600    AliAODJetEventBackground* externalBackground = 0;
 601    if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
 602       externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
 603       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 604    }
 605    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
 606      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
 607       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 608    }
 609
 610     if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
 611       externalBackground =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
 612       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 613     }
 614
 615    Float_t rho = 0;
 616
 617    if(fFlagRandom==0){
 618      if(externalBackground)rho = externalBackground->GetBackground(0);}
 619    if(fFlagRandom==1){
 620       if(externalBackground)rho = externalBackground->GetBackground(2);}
 621
 622    // fetch jets
 623    TClonesArray *aodJets[2];
 624    aodJets[0]=0;
 625    if(fAODOut&&!aodJets[0]){
 626    aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
 627    aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));  }
 628    if(fAODExtension && !aodJets[0]){
 629    aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
 630    aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
 631      if(fAODIn&&!aodJets[0]){
 632    aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
 633    aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data()));  }
 634
 635
 636    //Double_t ptsub[aodJets[0]->GetEntriesFast()];
 637    //Int_t inord[aodJets[0]->GetEntriesFast()];
 638    //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
 639    //  ptsub[n]=0;
 640    //  inord[n]=0;}
 641
 642    TList ParticleList;
 643    Int_t nT = GetListOfTracks(&ParticleList);
 644      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
 645       fListJets[iJetType]->Clear();
 646       if (!aodJets[iJetType]) continue;
 647       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
 648       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 649          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 650          if (jet) fListJets[iJetType]->Add(jet);
 651          // if(iJetType==0){
 652          // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
 653       }}
 654
 655    Double_t etabig=0;
 656    Double_t ptbig=0;
 657    Double_t areabig=0;
 658    Double_t phibig=0.;
 659    Double_t etasmall=0;
 660    Double_t ptsmall=0;
 661    Double_t areasmall=0;
 662    Double_t phismall=0.;
 663
 664
 665
 666    Int_t iCount=0;
 667    Int_t trigJet=-1;
 668    Int_t trigBBTrack=-1;
 669    Int_t trigInTrack=-1;
 670    fRPAngle = aod->GetHeader()->GetEventplane();
 671
 672    AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);
 673    if(!partback){
 674    PostData(1, fOutputList);
 675    return;}
 676
 677
 678    //for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
 679    //if(fFlagOnlyHardest!=0){if(tt!=nT) continue;}
 680    //AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);
 681
 682    fh2Ntriggers->Fill(centValue,partback->Pt());
 683    Double_t phiBinT = RelativePhi(partback->Phi(),fRPAngle);
 684    if(centValue<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
 685    if(centValue<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
 686    Double_t accep=2.*TMath::Pi()*1.8;
 687    Int_t injet4=0;
 688    Int_t injet=0;
 689
 690    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
 691            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
 692            etabig  = jetbig->Eta();
 693            phibig  = jetbig->Phi();
 694            ptbig   = jetbig->Pt();
 695            if(ptbig==0) continue;
 696            Double_t phiBin = RelativePhi(phibig,fRPAngle);
 697            areabig = jetbig->EffectiveAreaCharged();
 698            Double_t ptcorr=ptbig-rho*areabig;
 699            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
 700            if(areabig>=0.07) injet=injet+1;
 701            if(areabig>=0.4) injet4=injet4+1;
 702            Double_t dphi=RelativePhi(partback->Phi(),phibig);
 703
 704            if(fFlagEtaBkg==1){
 705            Double_t etadif= partback->Eta()-etabig;
 706            if(TMath::Abs(etadif)<=0.5){
 707
 708            if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
 709            if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}}
 710            if(fFlagEtaBkg==0){
 711            if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
 712            if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}
 713
 714
 715            if(fFlagJetHadron==0){
 716            if(fFlagPhiBkg!=0) if((TMath::Abs(dphi)<TMath::Pi()/2.-0.1)||(TMath::Abs(dphi)>TMath::Pi()/2.+0.1)) continue;
 717            if(fFlagPhiBkg==0) if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;}
 718
 719            if(fFlagJetHadron!=0) if(TMath::Abs(dphi)>0.4) continue;
 720
 721
 722            if(centValue<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
 723            if(centValue<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
 724                    Double_t dismin=100.;
 725                    Double_t ptmax=-10.;
 726                    Int_t index1=-1;
 727                    Int_t index2=-1;
 728
 729            if(centValue<10.)  fh3spectriggeredC10->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
 730            if(centValue<20.)  fh3spectriggeredC20->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
 731            if(centValue>30. && centValue<60.)  fh3spectriggeredC3060->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
 732
 733                    if(ptcorr<=0) continue;
 734
 735                        AliAODTrack* leadtrack=0;
 736                        Int_t ippt=0;
 737                        Double_t ppt=-10;
 738                        if(fFlagJetHadron==0){
 739                        TRefArray *genTrackList = jetbig->GetRefTracks();
 740                        Int_t nTracksGenJet = genTrackList->GetEntriesFast();
 741                        AliAODTrack* genTrack;
 742                        for(Int_t ir=0; ir<nTracksGenJet; ++ir){
 743                        genTrack = (AliAODTrack*)(genTrackList->At(ir));
 744                        if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
 745                        ippt=ir;}}
 746                        leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
 747                        if(!leadtrack) continue;}
 748
 749                        AliVParticle* leadtrackb=0;
 750                        if(fFlagJetHadron!=0){
 751                          Int_t nTb = GetHardestTrackBackToJet(jetbig);
 752                          leadtrackb = (AliVParticle*)ParticleList.At(nTb);
 753                          if(!leadtrackb) continue;
 754                        }
 755
 756
 757
 758
 759
 760                         //store one trigger info
 761                         if(iCount==0){
 762                           trigJet=i;
 763                           trigBBTrack=nT;
 764                           trigInTrack=ippt;
 765                           iCount=iCount+1;}
 766
 767
 768                   if(fCheckMethods){
 769                   for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
 770                   AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
 771                   etasmall  = jetsmall->Eta();
 772                   phismall = jetsmall->Phi();
 773                   ptsmall   = jetsmall->Pt();
 774                   areasmall = jetsmall->EffectiveAreaCharged();
 775                   Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
 776                   tmpDeltaR=TMath::Sqrt(tmpDeltaR);
 777                      //Fraction in the jet core
 778                     if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
 779                     index2=j;}
 780                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
 781                       index1=j;}} //en of loop over R=0.2 jets
 782                 //method1:most concentric jet=core
 783                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
 784                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 785                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 786                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 787                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
 788                 //method2:hardest contained jet=core
 789                 if(index2!=-1){
 790                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
 791                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 792                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 793                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 794                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}
 795                   if(centValue<10) fh2Ntriggers2C10->Fill(leadtrack->Pt(),partback->Pt());
 796                   if(centValue<20) fh2Ntriggers2C20->Fill(leadtrack->Pt(),partback->Pt());
 797          if(fDoEventMixing==0 && fFlagOnlyRecoil==0){
 798          for(int it = 0;it<ParticleList.GetEntries();++it){
 799           AliVParticle *part = (AliVParticle*)ParticleList.At(it);
 800           Double_t deltaR = jetbig->DeltaR(part);
 801           Double_t deltaEta = etabig-part->Eta();
 802
 803           Double_t deltaPhi=phibig-part->Phi();
 804           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
 805           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
 806           Double_t pTcont=0;
 807           if(fFlagJetHadron==0) pTcont=leadtrack->Pt();
 808           if(fFlagJetHadron!=0) pTcont=leadtrackb->Pt();
 809            Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,pTcont,partback->Pt()};
 810            fhnDeltaR->Fill(jetEntries);}
 811
 812
 813           }
 814
 815
 816           //end of track loop, we only do it if EM is switched off
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826    }
 827    if(injet>0) fh3JetDensity->Fill(ParticleList.GetEntries(),injet/accep,partback->Pt());
 828    if(injet4>0)fh3JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep,partback->Pt());
 829           //end of jet loop
 830
 831    //}
 832
 833
 834           if(fDoEventMixing>0){
 835             //check before if the trigger exists
 836             // fTrigBuffer[i][0] = zvtx
 837             // fTrigBuffer[i][1] = phi
 838             // fTrigBuffer[i][2] = eta
 839             // fTrigBuffer[i][3] = pt_jet
 840             // fTrigBuffer[i][4] = pt_trig
 841             // fTrigBuffer[i][5]= centrality
 842             if(fTindex==10) fTindex=0;
 843             if(fTrigBuffer[fTindex][3]>0){
 844             if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
 845             if (TMath::Abs(fTrigBuffer[fTindex][5]-centValue<5)){
 846
 847                         for(int it = 0;it<nT;++it){
 848                         AliVParticle *part = (AliVParticle*)ParticleList.At(it);
 849                         Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
 850                         Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
 851                         Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
 852                         if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
 853                         if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
 854                         Double_t triggerEntries[7] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4]};
 855                         fhnMixedEvents->Fill(triggerEntries);
 856                         }
 857                         fNevents=fNevents+1;
 858                         if(fNevents==10) fTindex=fTindex+1;
 859             }}}
 860
 861                if(fTindex==10&&fNevents==10) fCountAgain=0;
 862
 863                // Copy the triggers from the current event into the buffer.
 864                //again, only if the trigger exists:
 865                if(fCountAgain==0){
 866                 if(trigJet>-1){
 867                 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet));                      AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);
 868                 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
 869                 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
 870                 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
 871                 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
 872                 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
 873                 fTrigBuffer[fTrigBufferIndex][5] = centValue;
 874                 fTrigBufferIndex++;
 875                 if(fTrigBufferIndex==9) {fTrigBufferIndex=0;
 876                                          fCountAgain=1;}
 877                 }
 878                }
 879
 880           }
 881
 882
 883
 884      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
 885
 886      //tracks up to R=0.8 distant from the jet axis
 887    //   if(fAngStructCloseTracks==1){
 888    //    TList CloseTrackList;
 889    //    Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
 890    //    Double_t difR=0.04;
 891    //    for(Int_t l=0;l<15;l++){
 892    //    Double_t rr=l*0.1+0.1;
 893    //     for(int it = 0;it<nn;++it){
 894    //         AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
 895    //         for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){
 896    //         AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);
 897    //         Double_t ptm=part1->Pt();
 898    //         Double_t ptn=part2->Pt();
 899    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
 900    //                    Rnm=TMath::Sqrt(Rnm);
 901    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
 902    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
 903    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
 904    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
 905    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
 906    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
 907    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
 908    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
 909    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
 910    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
 911    //   }
 912
 913    //   //only jet constituents
 914    //    if(fAngStructCloseTracks==2){
 915
 916    //    Double_t difR=0.04;
 917    //    for(Int_t l=0;l<15;l++){
 918    //    Double_t rr=l*0.1+0.1;
 919
 920
 921    //    AliAODTrack* part1;
 922    //    AliAODTrack* part2;
 923
 924    //        TRefArray *genTrackListb = jetbig->GetRefTracks();
 925    //        Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
 926
 927
 928
 929    //        for(Int_t it=0; it<nTracksGenJetb; ++it){
 930    //           part1 = (AliAODTrack*)(genTrackListb->At(it));
 931    //         for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
 932    //           part2 = (AliAODTrack*)(genTrackListb->At(itu));
 933    //         Double_t ptm=part1->Pt();
 934    //         Double_t ptn=part2->Pt();
 935    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
 936    //                    Rnm=TMath::Sqrt(Rnm);
 937    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
 938    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
 939    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
 940    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
 941    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
 942    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
 943    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
 944    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
 945    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
 946    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
 947    // }
 948      // //end loop over R=0.4 jets
 949      // if(fAngStructCloseTracks>0){
 950      // for(Int_t l=0;l<15;l++){
 951      // Double_t rr=l*0.1+0.1;
 952      //    if(down1[l]!=0){
 953      //         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
 954      //    if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
 955      //    if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
 956      //    if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
 957      //    if(down2[l]!=0){
 958      //         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
 959      //    if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
 960      //    if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
 961      //    if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
 962      //    if(down3[l]!=0){
 963      //         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
 964      //    if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
 965      //    if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
 966      //    if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
 967      //    if(down4[l]!=0){
 968      //         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
 969      //    if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
 970      //    if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
 971      //    if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
 972
 973
 974
 975
 976
 977
 978
 979    PostData(1, fOutputList);
 980 }
 981
 982 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
 983 {
 984    // Draw result to the screen
 985    // Called once at the end of the query
 986
 987    if (!GetOutputData(1))
 988    return;
 989 }
 990
 991
 992
 993
 994
 995
 996 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
 997
 998       Int_t iCount = 0;
 999       AliAODEvent *aod = 0;
1000
1001      if(!fESD)aod = fAODIn;
1002      else aod = fAODOut;
1003       Int_t index=-1;
1004      Double_t ptmax=-10;
1005
1006
1007
1008      for(int it = 0;it < aod->GetNumberOfTracks();++it){
1009       AliAODTrack *tr = aod->GetTrack(it);
1010       Bool_t bGood = false;
1011       if(fFilterType == 0)bGood = true;
1012       else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1013       else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1014      if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1015       if(bGood==false) continue;
1016       if(TMath::Abs(tr->Eta())>0.9)continue;
1017       if(tr->Pt()<0.15)continue;
1018       list->Add(tr);
1019       iCount++;
1020       if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
1021         if(tr->TestFilterBit(fFilterMaskBestPt)){
1022           if(tr->Pt()>ptmax){
1023             ptmax=tr->Pt();
1024             index=iCount-1;
1025           }
1026         }
1027       }
1028       else{
1029         if(tr->Pt()>ptmax){
1030           ptmax=tr->Pt();
1031           index=iCount-1;
1032         }
1033       }
1034      }
1035
1036
1037     // else if (type == kTrackAODMCCharged) {
1038     // TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1039     // if(!tca)return iCount;
1040     // for(int it = 0;it < tca->GetEntriesFast();++it){
1041     //   AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1042     //   if(!part)continue;
1043     //   if(part->Pt()<0.15)continue;
1044     //   if(!part->IsPhysicalPrimary())continue;
1045     //   if(part->Charge()==0)continue;
1046     //   if(TMath::Abs(part->Eta())>0.9)continue;
1047     //   list->Add(part);
1048     //   iCount++;
1049     //   if(part->Pt()>ptmax){ ptmax=part->Pt();
1050     //  index=iCount-1;}}}
1051       return index;
1052
1053 }
1054
1055    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
1056
1057     AliAODEvent *aod = 0;
1058     if(!fESD)aod = fAODIn;
1059     else aod = fAODOut;
1060     Int_t index=-1;
1061     Double_t ptmax=-10;
1062     Double_t dphi=0;
1063     Double_t dif=0;
1064     Int_t iCount=0;
1065     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1066       AliAODTrack *tr = aod->GetTrack(it);
1067       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1068       if(TMath::Abs(tr->Eta())>0.9)continue;
1069       if(tr->Pt()<0.15)continue;
1070       iCount=iCount+1;
1071       dphi=RelativePhi(tr->Phi(),jetbig->Phi());
1072       if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
1073       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
1074       index=iCount-1;
1075       dif=dphi;  }}
1076
1077       return index;
1078
1079    }
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1090
1091     Int_t iCount = 0;
1092       AliAODEvent *aod = 0;
1093      if(!fESD)aod = fAODIn;
1094      else aod = fAODOut;
1095
1096       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1097       AliAODTrack *tr = aod->GetTrack(it);
1098       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1099       if(TMath::Abs(tr->Eta())>0.9)continue;
1100       if(tr->Pt()<0.15)continue;
1101       Double_t disR=jetbig->DeltaR(tr);
1102       if(disR>0.8)  continue;
1103       list->Add(tr);
1104       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1105       iCount++;
1106     }
1107
1108    list->Sort();
1109    return iCount;
1110
1111 }
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1124 {
1125
1126    Int_t nInputTracks = 0;
1127      AliAODEvent *aod = 0;
1128      if(!fESD)aod = fAODIn;
1129      else aod = fAODOut;
1130    TString jbname(fJetBranchName[1]);
1131    //needs complete event, use jets without background subtraction
1132    for(Int_t i=1; i<=3; ++i){
1133       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1134    }
1135    // use only HI event
1136    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1137    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1138
1139    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1140    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
1141    if(!tmpAODjets){
1142       Printf("Jet branch %s not found", jbname.Data());
1143       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1144       return -1;
1145    }
1146
1147    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1148       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1149       if(!jet) continue;
1150       TRefArray *trackList = jet->GetRefTracks();
1151       Int_t nTracks = trackList->GetEntriesFast();
1152       nInputTracks += nTracks;
1153       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1154    }
1155    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1156
1157    return nInputTracks;
1158 }
1159
1160
1161
1162 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1163
1164   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1165   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1166   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1167   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1168   double dphi = mphi-vphi;
1169   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1170   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1171   return dphi;//dphi in [-Pi, Pi]
1172 }
1173
1174 Int_t AliAnalysisTaskJetCore::GetPhiBin(Double_t phi)
1175 {
1176     Int_t phibin=-1;
1177     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1178     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1179     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1180     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1181     return phibin;
1182 }
1183
1184
1185
1186
1187 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1188 {
1189    // generate new THnSparseF, axes are defined in GetDimParams()
1190
1191    Int_t count = 0;
1192    UInt_t tmp = entries;
1193    while(tmp!=0){
1194       count++;
1195       tmp = tmp &~ -tmp;  // clear lowest bit
1196    }
1197
1198    TString hnTitle(name);
1199    const Int_t dim = count;
1200    Int_t nbins[dim];
1201    Double_t xmin[dim];
1202    Double_t xmax[dim];
1203
1204    Int_t i=0;
1205    Int_t c=0;
1206    while(c<dim && i<32){
1207       if(entries&(1<<i)){
1208
1209          TString label("");
1210          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1211          hnTitle += Form(";%s",label.Data());
1212          c++;
1213       }
1214
1215       i++;
1216    }
1217    hnTitle += ";";
1218
1219    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1220 }
1221
1222 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1223 {
1224    // stores label and binning of axis for THnSparse
1225
1226    const Double_t pi = TMath::Pi();
1227
1228    switch(iEntry){
1229
1230    case 0:
1231       label = "V0 centrality (%)";
1232
1233          nbins = 10;
1234          xmin = 0.;
1235          xmax = 100.;
1236          break;
1237
1238
1239    case 1:
1240       label = "corrected jet pt";
1241          nbins = 20;
1242          xmin = 0.;
1243          xmax = 200.;
1244           break;
1245
1246
1247    case 2:
1248       label = "track pT";
1249
1250          nbins = 9;
1251          xmin = 0.;
1252          xmax = 150;
1253          break;
1254
1255
1256     case 3:
1257       label = "deltaR";
1258       nbins = 15;
1259       xmin = 0.;
1260       xmax = 1.5;
1261       break;
1262
1263
1264
1265    case 4:
1266       label = "deltaEta";
1267       nbins = 8;
1268       xmin = -1.6;
1269       xmax = 1.6;
1270       break;
1271
1272
1273   case 5:
1274       label = "deltaPhi";
1275       nbins = 90;
1276       xmin = -0.5*pi;
1277       xmax = 1.5*pi;
1278       break;
1279
1280
1281
1282     case 6:
1283       label = "leading track";
1284       nbins = 13;
1285       xmin = 0;
1286       xmax = 50;
1287       break;
1288
1289      case 7:
1290
1291       label = "trigger track";
1292       nbins =10;
1293       xmin = 0;
1294       xmax = 50;
1295       break;
1296
1297
1298
1299
1300
1301
1302
1303
1304    }
1305
1306 }
1307