| | 1 | // ************************************** |
| | 2 | // Task used for the correction of determiantion of reconstructed jet spectra |
| | 3 | // Compares input (gen) and output (rec) jets |
| | 4 | // ******************************************* |
| | 5 | |
| | 6 | |
| | 7 | /************************************************************************** |
| | 8 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * |
| | 9 | * * |
| | 10 | * Author: The ALICE Off-line Project. * |
| | 11 | * Contributors are mentioned in the code where appropriate. * |
| | 12 | * * |
| | 13 | * Permission to use, copy, modify and distribute this software and its * |
| | 14 | * documentation strictly for non-commercial purposes is hereby granted * |
| | 15 | * without fee, provided that the above copyright notice appears in all * |
| | 16 | * copies and that both the copyright notice and this permission notice * |
| | 17 | * appear in the supporting documentation. The authors make no claims * |
| | 18 | * about the suitability of this software for any purpose. It is * |
| | 19 | * provided "as is" without express or implied warranty. * |
| | 20 | **************************************************************************/ |
| | 21 | |
| | 22 | |
| | 23 | #include <TROOT.h> |
| | 24 | #include <TRandom3.h> |
| | 25 | #include <TSystem.h> |
| | 26 | #include <TInterpreter.h> |
| | 27 | #include <TChain.h> |
| | 28 | #include <TRefArray.h> |
| | 29 | #include <TFile.h> |
| | 30 | #include <TKey.h> |
| | 31 | #include <TH1F.h> |
| | 32 | #include <TH2F.h> |
| | 33 | #include <TH3F.h> |
| | 34 | #include <TProfile.h> |
| | 35 | #include <TF1.h> |
| | 36 | #include <TList.h> |
| | 37 | #include <TLorentzVector.h> |
| | 38 | #include <TClonesArray.h> |
| | 39 | #include "TDatabasePDG.h" |
| | 40 | |
| | 41 | #include "AliAnalysisTaskJetCluster.h" |
| | 42 | #include "AliAnalysisManager.h" |
| | 43 | #include "AliJetFinder.h" |
| | 44 | #include "AliJetHeader.h" |
| | 45 | #include "AliJetReader.h" |
| | 46 | #include "AliESDEvent.h" |
| | 47 | #include "AliAODEvent.h" |
| | 48 | #include "AliAODHandler.h" |
| | 49 | #include "AliAODExtension.h" |
| | 50 | #include "AliAODTrack.h" |
| | 51 | #include "AliAODJet.h" |
| | 52 | #include "AliAODMCParticle.h" |
| | 53 | #include "AliMCEventHandler.h" |
| | 54 | #include "AliMCEvent.h" |
| | 55 | #include "AliStack.h" |
| | 56 | #include "AliGenPythiaEventHeader.h" |
| | 57 | #include "AliJetKineReaderHeader.h" |
| | 58 | #include "AliGenCocktailEventHeader.h" |
| | 59 | #include "AliInputEventHandler.h" |
| | 60 | #include "AliAODJetEventBackground.h" |
| | 61 | |
| | 62 | #include "fastjet/PseudoJet.hh" |
| | 63 | #include "fastjet/ClusterSequenceArea.hh" |
| | 64 | #include "fastjet/AreaDefinition.hh" |
| | 65 | #include "fastjet/JetDefinition.hh" |
| | 66 | // get info on how fastjet was configured |
| | 67 | #include "fastjet/config.h" |
| | 68 | |
| | 69 | using std::vector; |
| | 70 | |
| | 71 | ClassImp(AliAnalysisTaskJetCluster) |
| | 72 | |
| | 73 | AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){ |
| | 74 | // |
| | 75 | // Destructor |
| | 76 | // |
| | 77 | |
| | 78 | delete fRef; |
| | 79 | delete fRandom; |
| | 80 | |
| | 81 | if(fTCAJetsOut)fTCAJetsOut->Delete(); |
| | 82 | delete fTCAJetsOut; |
| | 83 | |
| | 84 | if(fTCAJetsOutRan)fTCAJetsOutRan->Delete(); |
| | 85 | delete fTCAJetsOutRan; |
| | 86 | |
| | 87 | if(fTCARandomConesOut)fTCARandomConesOut->Delete(); |
| | 88 | delete fTCARandomConesOut; |
| | 89 | |
| | 90 | if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete(); |
| | 91 | delete fTCARandomConesOutRan; |
| | 92 | |
| | 93 | if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset(); |
| | 94 | delete fAODJetBackgroundOut; |
| | 95 | } |
| | 96 | |
| | 97 | AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): |
| | 98 | AliAnalysisTaskSE(), |
| | 99 | fAOD(0x0), |
| | 100 | fAODExtension(0x0), |
| | 101 | fRef(new TRefArray), |
| | 102 | fUseAODTrackInput(kFALSE), |
| | 103 | fUseAODMCInput(kFALSE), |
| | 104 | fUseBackgroundCalc(kFALSE), |
| | 105 | fEventSelection(kFALSE), |
| | 106 | fFilterMask(0), |
| | 107 | fFilterMaskBestPt(0), |
| | 108 | fFilterType(0), |
| | 109 | fJetTypes(kJet), |
| | 110 | fTrackTypeRec(kTrackUndef), |
| | 111 | fTrackTypeGen(kTrackUndef), |
| | 112 | fNSkipLeadingRan(0), |
| | 113 | fNSkipLeadingCone(0), |
| | 114 | fNRandomCones(0), |
| | 115 | fAvgTrials(1), |
| | 116 | fExternalWeight(1), |
| | 117 | fTrackEtaWindow(0.9), |
| | 118 | fRecEtaWindow(0.5), |
| | 119 | fTrackPtCut(0.), |
| | 120 | fJetOutputMinPt(0.150), |
| | 121 | fMaxTrackPtInJet(100.), |
| | 122 | fJetTriggerPtCut(0), |
| | 123 | fVtxZCut(8), |
| | 124 | fVtxR2Cut(1), |
| | 125 | fCentCutUp(0), |
| | 126 | fCentCutLo(0), |
| | 127 | fNonStdBranch(""), |
| | 128 | fBackgroundBranch(""), |
| | 129 | fNonStdFile(""), |
| | 130 | fMomResH1(0x0), |
| | 131 | fMomResH2(0x0), |
| | 132 | fMomResH3(0x0), |
| | 133 | fMomResH1Fit(0x0), |
| | 134 | fMomResH2Fit(0x0), |
| | 135 | fMomResH3Fit(0x0), |
| | 136 | fhEffH1(0x0), |
| | 137 | fhEffH2(0x0), |
| | 138 | fhEffH3(0x0), |
| | 139 | fUseTrPtResolutionSmearing(kFALSE), |
| | 140 | fUseDiceEfficiency(kFALSE), |
| | 141 | fDiceEfficiencyMinPt(0.), |
| | 142 | fUseTrPtResolutionFromOADB(kFALSE), |
| | 143 | fUseTrEfficiencyFromOADB(kFALSE), |
| | 144 | fPathTrPtResolution(""), |
| | 145 | fPathTrEfficiency(""), |
| | 146 | fChangeEfficiencyFraction(0.), |
| | 147 | fRparam(1.0), |
| | 148 | fAlgorithm(fastjet::kt_algorithm), |
| | 149 | fStrategy(fastjet::Best), |
| | 150 | fRecombScheme(fastjet::BIpt_scheme), |
| | 151 | fAreaType(fastjet::active_area), |
| | 152 | fGhostArea(0.01), |
| | 153 | fActiveAreaRepeats(1), |
| | 154 | fGhostEtamax(1.5), |
| | 155 | fTCAJetsOut(0x0), |
| | 156 | fTCAJetsOutRan(0x0), |
| | 157 | fTCARandomConesOut(0x0), |
| | 158 | fTCARandomConesOutRan(0x0), |
| | 159 | fAODJetBackgroundOut(0x0), |
| | 160 | fRandom(0), |
| | 161 | fh1Xsec(0x0), |
| | 162 | fh1Trials(0x0), |
| | 163 | fh1PtHard(0x0), |
| | 164 | fh1PtHardNoW(0x0), |
| | 165 | fh1PtHardTrials(0x0), |
| | 166 | fh1NJetsRec(0x0), |
| | 167 | fh1NConstRec(0x0), |
| | 168 | fh1NConstLeadingRec(0x0), |
| | 169 | fh1PtJetsRecIn(0x0), |
| | 170 | fh1PtJetsLeadingRecIn(0x0), |
| | 171 | fh1PtJetConstRec(0x0), |
| | 172 | fh1PtJetConstLeadingRec(0x0), |
| | 173 | fh1PtTracksRecIn(0x0), |
| | 174 | fh1PtTracksLeadingRecIn(0x0), |
| | 175 | fh1NJetsRecRan(0x0), |
| | 176 | fh1NConstRecRan(0x0), |
| | 177 | fh1PtJetsLeadingRecInRan(0x0), |
| | 178 | fh1NConstLeadingRecRan(0x0), |
| | 179 | fh1PtJetsRecInRan(0x0), |
| | 180 | fh1PtTracksGenIn(0x0), |
| | 181 | fh1Nch(0x0), |
| | 182 | fh1CentralityPhySel(0x0), |
| | 183 | fh1Centrality(0x0), |
| | 184 | fh1CentralitySelect(0x0), |
| | 185 | fh1ZPhySel(0x0), |
| | 186 | fh1Z(0x0), |
| | 187 | fh1ZSelect(0x0), |
| | 188 | fh2NRecJetsPt(0x0), |
| | 189 | fh2NRecTracksPt(0x0), |
| | 190 | fh2NConstPt(0x0), |
| | 191 | fh2NConstLeadingPt(0x0), |
| | 192 | fh2JetPhiEta(0x0), |
| | 193 | fh2LeadingJetPhiEta(0x0), |
| | 194 | fh2JetEtaPt(0x0), |
| | 195 | fh2LeadingJetEtaPt(0x0), |
| | 196 | fh2TrackEtaPt(0x0), |
| | 197 | fh2LeadingTrackEtaPt(0x0), |
| | 198 | fh2JetsLeadingPhiEta(0x0), |
| | 199 | fh2JetsLeadingPhiPt(0x0), |
| | 200 | fh2TracksLeadingPhiEta(0x0), |
| | 201 | fh2TracksLeadingPhiPt(0x0), |
| | 202 | fh2TracksLeadingJetPhiPt(0x0), |
| | 203 | fh2JetsLeadingPhiPtW(0x0), |
| | 204 | fh2TracksLeadingPhiPtW(0x0), |
| | 205 | fh2TracksLeadingJetPhiPtW(0x0), |
| | 206 | fh2NRecJetsPtRan(0x0), |
| | 207 | fh2NConstPtRan(0x0), |
| | 208 | fh2NConstLeadingPtRan(0x0), |
| | 209 | fh2PtNch(0x0), |
| | 210 | fh2PtNchRan(0x0), |
| | 211 | fh2PtNchN(0x0), |
| | 212 | fh2PtNchNRan(0x0), |
| | 213 | fh2TracksLeadingJetPhiPtRan(0x0), |
| | 214 | fh2TracksLeadingJetPhiPtWRan(0x0), |
| | 215 | fh2PtGenPtSmeared(0x0), |
| | 216 | fp1Efficiency(0x0), |
| | 217 | fp1PtResolution(0x0), |
| | 218 | fHistList(0x0) |
| | 219 | { |
| | 220 | // |
| | 221 | // Constructor |
| | 222 | // |
| | 223 | |
| | 224 | for(int i = 0;i<3;i++){ |
| | 225 | fh1BiARandomCones[i] = 0; |
| | 226 | fh1BiARandomConesRan[i] = 0; |
| | 227 | } |
| | 228 | for(int i = 0;i<kMaxCent;i++){ |
| | 229 | fh2JetsLeadingPhiPtC[i] = 0; |
| | 230 | fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet |
| | 231 | fh2TracksLeadingJetPhiPtC[i] = 0; |
| | 232 | fh2TracksLeadingJetPhiPtWC[i] = 0; |
| | 233 | } |
| | 234 | } |
| | 235 | |
| | 236 | AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): |
| | 237 | AliAnalysisTaskSE(name), |
| | 238 | fAOD(0x0), |
| | 239 | fAODExtension(0x0), |
| | 240 | fRef(new TRefArray), |
| | 241 | fUseAODTrackInput(kFALSE), |
| | 242 | fUseAODMCInput(kFALSE), |
| | 243 | fUseBackgroundCalc(kFALSE), |
| | 244 | fEventSelection(kFALSE), fFilterMask(0), |
| | 245 | fFilterMaskBestPt(0), |
| | 246 | fFilterType(0), |
| | 247 | fJetTypes(kJet), |
| | 248 | fTrackTypeRec(kTrackUndef), |
| | 249 | fTrackTypeGen(kTrackUndef), |
| | 250 | fNSkipLeadingRan(0), |
| | 251 | fNSkipLeadingCone(0), |
| | 252 | fNRandomCones(0), |
| | 253 | fAvgTrials(1), |
| | 254 | fExternalWeight(1), |
| | 255 | fTrackEtaWindow(0.9), |
| | 256 | fRecEtaWindow(0.5), |
| | 257 | fTrackPtCut(0.), |
| | 258 | fJetOutputMinPt(0.150), |
| | 259 | fMaxTrackPtInJet(100.), |
| | 260 | fJetTriggerPtCut(0), |
| | 261 | fVtxZCut(8), |
| | 262 | fVtxR2Cut(1), |
| | 263 | fCentCutUp(0), |
| | 264 | fCentCutLo(0), |
| | 265 | fNonStdBranch(""), |
| | 266 | fBackgroundBranch(""), |
| | 267 | fNonStdFile(""), |
| | 268 | fMomResH1(0x0), |
| | 269 | fMomResH2(0x0), |
| | 270 | fMomResH3(0x0), |
| | 271 | fMomResH1Fit(0x0), |
| | 272 | fMomResH2Fit(0x0), |
| | 273 | fMomResH3Fit(0x0), |
| | 274 | fhEffH1(0x0), |
| | 275 | fhEffH2(0x0), |
| | 276 | fhEffH3(0x0), |
| | 277 | fUseTrPtResolutionSmearing(kFALSE), |
| | 278 | fUseDiceEfficiency(kFALSE), |
| | 279 | fDiceEfficiencyMinPt(0.), |
| | 280 | fUseTrPtResolutionFromOADB(kFALSE), |
| | 281 | fUseTrEfficiencyFromOADB(kFALSE), |
| | 282 | fPathTrPtResolution(""), |
| | 283 | fPathTrEfficiency(""), |
| | 284 | fChangeEfficiencyFraction(0.), |
| | 285 | fRparam(1.0), |
| | 286 | fAlgorithm(fastjet::kt_algorithm), |
| | 287 | fStrategy(fastjet::Best), |
| | 288 | fRecombScheme(fastjet::BIpt_scheme), |
| | 289 | fAreaType(fastjet::active_area), |
| | 290 | fGhostArea(0.01), |
| | 291 | fActiveAreaRepeats(1), |
| | 292 | fGhostEtamax(1.5), |
| | 293 | fTCAJetsOut(0x0), |
| | 294 | fTCAJetsOutRan(0x0), |
| | 295 | fTCARandomConesOut(0x0), |
| | 296 | fTCARandomConesOutRan(0x0), |
| | 297 | fAODJetBackgroundOut(0x0), |
| | 298 | fRandom(0), |
| | 299 | fh1Xsec(0x0), |
| | 300 | fh1Trials(0x0), |
| | 301 | fh1PtHard(0x0), |
| | 302 | fh1PtHardNoW(0x0), |
| | 303 | fh1PtHardTrials(0x0), |
| | 304 | fh1NJetsRec(0x0), |
| | 305 | fh1NConstRec(0x0), |
| | 306 | fh1NConstLeadingRec(0x0), |
| | 307 | fh1PtJetsRecIn(0x0), |
| | 308 | fh1PtJetsLeadingRecIn(0x0), |
| | 309 | fh1PtJetConstRec(0x0), |
| | 310 | fh1PtJetConstLeadingRec(0x0), |
| | 311 | fh1PtTracksRecIn(0x0), |
| | 312 | fh1PtTracksLeadingRecIn(0x0), |
| | 313 | fh1NJetsRecRan(0x0), |
| | 314 | fh1NConstRecRan(0x0), |
| | 315 | fh1PtJetsLeadingRecInRan(0x0), |
| | 316 | fh1NConstLeadingRecRan(0x0), |
| | 317 | fh1PtJetsRecInRan(0x0), |
| | 318 | fh1PtTracksGenIn(0x0), |
| | 319 | fh1Nch(0x0), |
| | 320 | fh1CentralityPhySel(0x0), |
| | 321 | fh1Centrality(0x0), |
| | 322 | fh1CentralitySelect(0x0), |
| | 323 | fh1ZPhySel(0x0), |
| | 324 | fh1Z(0x0), |
| | 325 | fh1ZSelect(0x0), |
| | 326 | fh2NRecJetsPt(0x0), |
| | 327 | fh2NRecTracksPt(0x0), |
| | 328 | fh2NConstPt(0x0), |
| | 329 | fh2NConstLeadingPt(0x0), |
| | 330 | fh2JetPhiEta(0x0), |
| | 331 | fh2LeadingJetPhiEta(0x0), |
| | 332 | fh2JetEtaPt(0x0), |
| | 333 | fh2LeadingJetEtaPt(0x0), |
| | 334 | fh2TrackEtaPt(0x0), |
| | 335 | fh2LeadingTrackEtaPt(0x0), |
| | 336 | fh2JetsLeadingPhiEta(0x0), |
| | 337 | fh2JetsLeadingPhiPt(0x0), |
| | 338 | fh2TracksLeadingPhiEta(0x0), |
| | 339 | fh2TracksLeadingPhiPt(0x0), |
| | 340 | fh2TracksLeadingJetPhiPt(0x0), |
| | 341 | fh2JetsLeadingPhiPtW(0x0), |
| | 342 | fh2TracksLeadingPhiPtW(0x0), |
| | 343 | fh2TracksLeadingJetPhiPtW(0x0), |
| | 344 | fh2NRecJetsPtRan(0x0), |
| | 345 | fh2NConstPtRan(0x0), |
| | 346 | fh2NConstLeadingPtRan(0x0), |
| | 347 | fh2PtNch(0x0), |
| | 348 | fh2PtNchRan(0x0), |
| | 349 | fh2PtNchN(0x0), |
| | 350 | fh2PtNchNRan(0x0), |
| | 351 | fh2TracksLeadingJetPhiPtRan(0x0), |
| | 352 | fh2TracksLeadingJetPhiPtWRan(0x0), |
| | 353 | fh2PtGenPtSmeared(0x0), |
| | 354 | fp1Efficiency(0x0), |
| | 355 | fp1PtResolution(0x0), |
| | 356 | fHistList(0x0) |
| | 357 | { |
| | 358 | // |
| | 359 | // named ctor |
| | 360 | // |
| | 361 | |
| | 362 | for(int i = 0;i<3;i++){ |
| | 363 | fh1BiARandomCones[i] = 0; |
| | 364 | fh1BiARandomConesRan[i] = 0; |
| | 365 | } |
| | 366 | for(int i = 0;i<kMaxCent;i++){ |
| | 367 | fh2JetsLeadingPhiPtC[i] = 0; |
| | 368 | fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet |
| | 369 | fh2TracksLeadingJetPhiPtC[i] = 0; |
| | 370 | fh2TracksLeadingJetPhiPtWC[i] = 0; |
| | 371 | } |
| | 372 | DefineOutput(1, TList::Class()); |
| | 373 | } |
| | 374 | |
| | 375 | |
| | 376 | |
| | 377 | Bool_t AliAnalysisTaskJetCluster::Notify() |
| | 378 | { |
| | 379 | // |
| | 380 | // Implemented Notify() to read the cross sections |
| | 381 | // and number of trials from pyxsec.root |
| | 382 | // |
| | 383 | return kTRUE; |
| | 384 | } |
| | 385 | |
| | 386 | void AliAnalysisTaskJetCluster::UserCreateOutputObjects() |
| | 387 | { |
| | 388 | |
| | 389 | // |
| | 390 | // Create the output container |
| | 391 | // |
| | 392 | |
| | 393 | fRandom = new TRandom3(0); |
| | 394 | |
| | 395 | |
| | 396 | // Connect the AOD |
| | 397 | |
| | 398 | |
| | 399 | if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n"); |
| | 400 | |
| | 401 | |
| | 402 | |
| | 403 | if(fNonStdBranch.Length()!=0) |
| | 404 | { |
| | 405 | // only create the output branch if we have a name |
| | 406 | // Create a new branch for jets... |
| | 407 | // -> cleared in the UserExec.... |
| | 408 | // here we can also have the case that the brnaches are written to a separate file |
| | 409 | |
| | 410 | if(fJetTypes&kJet){ |
| | 411 | fTCAJetsOut = new TClonesArray("AliAODJet", 0); |
| | 412 | fTCAJetsOut->SetName(fNonStdBranch.Data()); |
| | 413 | AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data()); |
| | 414 | } |
| | 415 | |
| | 416 | if(fJetTypes&kJetRan){ |
| | 417 | fTCAJetsOutRan = new TClonesArray("AliAODJet", 0); |
| | 418 | fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random")); |
| | 419 | if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) |
| | 420 | fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.))); |
| | 421 | AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data()); |
| | 422 | } |
| | 423 | |
| | 424 | if(fUseBackgroundCalc){ |
| | 425 | if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){ |
| | 426 | fAODJetBackgroundOut = new AliAODJetEventBackground(); |
| | 427 | fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())); |
| | 428 | if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) |
| | 429 | fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.))); |
| | 430 | |
| | 431 | AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data()); |
| | 432 | } |
| | 433 | } |
| | 434 | // create the branch for the random cones with the same R |
| | 435 | TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone); |
| | 436 | if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) |
| | 437 | cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone); |
| | 438 | |
| | 439 | if(fNRandomCones>0){ |
| | 440 | if(fJetTypes&kRC){ |
| | 441 | if(!AODEvent()->FindListObject(cName.Data())){ |
| | 442 | fTCARandomConesOut = new TClonesArray("AliAODJet", 0); |
| | 443 | fTCARandomConesOut->SetName(cName.Data()); |
| | 444 | AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data()); |
| | 445 | } |
| | 446 | } |
| | 447 | // create the branch with the random for the random cones on the random event |
| | 448 | if(fJetTypes&kRCRan){ |
| | 449 | cName = Form("%sRandomCone_random",fNonStdBranch.Data()); |
| | 450 | if(!AODEvent()->FindListObject(cName.Data())){ |
| | 451 | fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0); |
| | 452 | fTCARandomConesOutRan->SetName(cName.Data()); |
| | 453 | AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data()); |
| | 454 | } |
| | 455 | } |
| | 456 | } |
| | 457 | |
| | 458 | if(fNonStdFile.Length()!=0){ |
| | 459 | // |
| | 460 | // case that we have an AOD extension we need to fetch the jets from the extended output |
| | 461 | // we identify the extension aod event by looking for the branchname |
| | 462 | AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); |
| | 463 | // case that we have an AOD extension we need can fetch the background maybe from the extended output |
| | 464 | fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0); |
| | 465 | } |
| | 466 | } |
| | 467 | |
| | 468 | |
| | 469 | if(!fHistList)fHistList = new TList(); |
| | 470 | fHistList->SetOwner(); |
| | 471 | PostData(1, fHistList); // post data in any case once |
| | 472 | |
| | 473 | Bool_t oldStatus = TH1::AddDirectoryStatus(); |
| | 474 | TH1::AddDirectory(kFALSE); |
| | 475 | |
| | 476 | // |
| | 477 | // Histogram |
| | 478 | |
| | 479 | const Int_t nBinPt = 100; |
| | 480 | Double_t binLimitsPt[nBinPt+1]; |
| | 481 | for(Int_t iPt = 0;iPt <= nBinPt;iPt++){ |
| | 482 | if(iPt == 0){ |
| | 483 | binLimitsPt[iPt] = 0.0; |
| | 484 | } |
| | 485 | else {// 1.0 |
| | 486 | binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0; |
| | 487 | } |
| | 488 | } |
| | 489 | |
| | 490 | const Int_t nBinPhi = 90; |
| | 491 | Double_t binLimitsPhi[nBinPhi+1]; |
| | 492 | for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){ |
| | 493 | if(iPhi==0){ |
| | 494 | binLimitsPhi[iPhi] = -1.*TMath::Pi(); |
| | 495 | } |
| | 496 | else{ |
| | 497 | binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2; |
| | 498 | } |
| | 499 | } |
| | 500 | |
| | 501 | |
| | 502 | |
| | 503 | const Int_t nBinEta = 40; |
| | 504 | Double_t binLimitsEta[nBinEta+1]; |
| | 505 | for(Int_t iEta = 0;iEta<=nBinEta;iEta++){ |
| | 506 | if(iEta==0){ |
| | 507 | binLimitsEta[iEta] = -2.0; |
| | 508 | } |
| | 509 | else{ |
| | 510 | binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1; |
| | 511 | } |
| | 512 | } |
| | 513 | |
| | 514 | const int nChMax = 5000; |
| | 515 | |
| | 516 | fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1); |
| | 517 | fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>"); |
| | 518 | |
| | 519 | fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1); |
| | 520 | fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}"); |
| | 521 | |
| | 522 | |
| | 523 | fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5); |
| | 524 | fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5); |
| | 525 | |
| | 526 | fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5); |
| | 527 | fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5); |
| | 528 | fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5); |
| | 529 | fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5); |
| | 530 | |
| | 531 | |
| | 532 | fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt); |
| | 533 | fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt); |
| | 534 | fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt); |
| | 535 | |
| | 536 | fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 537 | fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 538 | fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 539 | fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 540 | fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 541 | fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 542 | fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt); |
| | 543 | fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt); |
| | 544 | fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt); |
| | 545 | fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5); |
| | 546 | |
| | 547 | fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5); |
| | 548 | fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5); |
| | 549 | fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5); |
| | 550 | |
| | 551 | fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25); |
| | 552 | fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25); |
| | 553 | fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25); |
| | 554 | |
| | 555 | fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 556 | fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 557 | fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 558 | // |
| | 559 | |
| | 560 | |
| | 561 | fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 562 | fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 563 | fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 564 | fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 565 | |
| | 566 | fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt); |
| | 567 | fh2PtNchRan = new TH2F("fh2PtNchRan","p_T of cluster vs. multiplicity ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt); |
| | 568 | fh2PtNchN = new TH2F("fh2PtNchN","p_T of cluster vs. multiplicity N weighted; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt); |
| | 569 | fh2PtNchNRan = new TH2F("fh2PtNchNRan","p_T of cluster vs. multiplicity N weighted ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt); |
| | 570 | |
| | 571 | |
| | 572 | |
| | 573 | fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta", |
| | 574 | nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta); |
| | 575 | fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta", |
| | 576 | nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta); |
| | 577 | |
| | 578 | fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}", |
| | 579 | nBinEta,binLimitsEta,nBinPt,binLimitsPt); |
| | 580 | fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}", |
| | 581 | nBinEta,binLimitsEta,nBinPt,binLimitsPt); |
| | 582 | |
| | 583 | fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}", |
| | 584 | nBinEta,binLimitsEta,nBinPt,binLimitsPt); |
| | 585 | fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}", |
| | 586 | nBinEta,binLimitsEta,nBinPt,binLimitsPt); |
| | 587 | |
| | 588 | |
| | 589 | |
| | 590 | fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta", |
| | 591 | nBinPhi,binLimitsPhi,nBinEta,binLimitsEta); |
| | 592 | fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 593 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 594 | fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta", |
| | 595 | nBinPhi,binLimitsPhi,nBinEta,binLimitsEta); |
| | 596 | fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 597 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 598 | fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 599 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 600 | fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 601 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 602 | |
| | 603 | fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 604 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 605 | fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)", |
| | 606 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 607 | |
| | 608 | fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 609 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 610 | fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 611 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 612 | |
| | 613 | //Detector level effects histos |
| | 614 | fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt); |
| | 615 | |
| | 616 | fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt); |
| | 617 | fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt); |
| | 618 | |
| | 619 | fHistList->Add(fh2PtGenPtSmeared); |
| | 620 | fHistList->Add(fp1Efficiency); |
| | 621 | fHistList->Add(fp1PtResolution); |
| | 622 | |
| | 623 | if(fNRandomCones>0&&fUseBackgroundCalc){ |
| | 624 | for(int i = 0;i<3;i++){ |
| | 625 | fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100); |
| | 626 | fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100); |
| | 627 | } |
| | 628 | } |
| | 629 | |
| | 630 | for(int i = 0;i < kMaxCent;i++){ |
| | 631 | fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1)); |
| | 632 | fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1)); |
| | 633 | fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1)); |
| | 634 | fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1)); |
| | 635 | } |
| | 636 | |
| | 637 | const Int_t saveLevel = 3; // large save level more histos |
| | 638 | if(saveLevel>0){ |
| | 639 | fHistList->Add(fh1Xsec); |
| | 640 | fHistList->Add(fh1Trials); |
| | 641 | |
| | 642 | fHistList->Add(fh1NJetsRec); |
| | 643 | fHistList->Add(fh1NConstRec); |
| | 644 | fHistList->Add(fh1NConstLeadingRec); |
| | 645 | fHistList->Add(fh1PtJetsRecIn); |
| | 646 | fHistList->Add(fh1PtJetsLeadingRecIn); |
| | 647 | fHistList->Add(fh1PtTracksRecIn); |
| | 648 | fHistList->Add(fh1PtTracksLeadingRecIn); |
| | 649 | fHistList->Add(fh1PtJetConstRec); |
| | 650 | fHistList->Add(fh1PtJetConstLeadingRec); |
| | 651 | fHistList->Add(fh1NJetsRecRan); |
| | 652 | fHistList->Add(fh1NConstRecRan); |
| | 653 | fHistList->Add(fh1PtJetsLeadingRecInRan); |
| | 654 | fHistList->Add(fh1NConstLeadingRecRan); |
| | 655 | fHistList->Add(fh1PtJetsRecInRan); |
| | 656 | fHistList->Add(fh1Nch); |
| | 657 | fHistList->Add(fh1Centrality); |
| | 658 | fHistList->Add(fh1CentralitySelect); |
| | 659 | fHistList->Add(fh1CentralityPhySel); |
| | 660 | fHistList->Add(fh1Z); |
| | 661 | fHistList->Add(fh1ZSelect); |
| | 662 | fHistList->Add(fh1ZPhySel); |
| | 663 | if(fNRandomCones>0&&fUseBackgroundCalc){ |
| | 664 | for(int i = 0;i<3;i++){ |
| | 665 | fHistList->Add(fh1BiARandomCones[i]); |
| | 666 | fHistList->Add(fh1BiARandomConesRan[i]); |
| | 667 | } |
| | 668 | } |
| | 669 | for(int i = 0;i < kMaxCent;i++){ |
| | 670 | fHistList->Add(fh2JetsLeadingPhiPtC[i]); |
| | 671 | fHistList->Add(fh2JetsLeadingPhiPtWC[i]); |
| | 672 | fHistList->Add(fh2TracksLeadingJetPhiPtC[i]); |
| | 673 | fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]); |
| | 674 | } |
| | 675 | |
| | 676 | fHistList->Add(fh2NRecJetsPt); |
| | 677 | fHistList->Add(fh2NRecTracksPt); |
| | 678 | fHistList->Add(fh2NConstPt); |
| | 679 | fHistList->Add(fh2NConstLeadingPt); |
| | 680 | fHistList->Add(fh2PtNch); |
| | 681 | fHistList->Add(fh2PtNchRan); |
| | 682 | fHistList->Add(fh2PtNchN); |
| | 683 | fHistList->Add(fh2PtNchNRan); |
| | 684 | fHistList->Add(fh2JetPhiEta); |
| | 685 | fHistList->Add(fh2LeadingJetPhiEta); |
| | 686 | fHistList->Add(fh2JetEtaPt); |
| | 687 | fHistList->Add(fh2LeadingJetEtaPt); |
| | 688 | fHistList->Add(fh2TrackEtaPt); |
| | 689 | fHistList->Add(fh2LeadingTrackEtaPt); |
| | 690 | fHistList->Add(fh2JetsLeadingPhiEta ); |
| | 691 | fHistList->Add(fh2JetsLeadingPhiPt); |
| | 692 | fHistList->Add(fh2TracksLeadingPhiEta); |
| | 693 | fHistList->Add(fh2TracksLeadingPhiPt); |
| | 694 | fHistList->Add(fh2TracksLeadingJetPhiPt); |
| | 695 | fHistList->Add(fh2JetsLeadingPhiPtW); |
| | 696 | fHistList->Add(fh2TracksLeadingPhiPtW); |
| | 697 | fHistList->Add(fh2TracksLeadingJetPhiPtW); |
| | 698 | fHistList->Add(fh2NRecJetsPtRan); |
| | 699 | fHistList->Add(fh2NConstPtRan); |
| | 700 | fHistList->Add(fh2NConstLeadingPtRan); |
| | 701 | fHistList->Add(fh2TracksLeadingJetPhiPtRan); |
| | 702 | fHistList->Add(fh2TracksLeadingJetPhiPtWRan); |
| | 703 | } |
| | 704 | |
| | 705 | // =========== Switch on Sumw2 for all histos =========== |
| | 706 | for (Int_t i=0; i<fHistList->GetEntries(); ++i) { |
| | 707 | TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i)); |
| | 708 | if (h1){ |
| | 709 | h1->Sumw2(); |
| | 710 | continue; |
| | 711 | } |
| | 712 | THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i)); |
| | 713 | if(hn)hn->Sumw2(); |
| | 714 | } |
| | 715 | TH1::AddDirectory(oldStatus); |
| | 716 | } |
| | 717 | |
| | 718 | void AliAnalysisTaskJetCluster::LocalInit() |
| | 719 | { |
| | 720 | // |
| | 721 | // Initialization |
| | 722 | // |
| | 723 | |
| | 724 | if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n"); |
| | 725 | |
| | 726 | if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB(); |
| | 727 | if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB(); |
| | 728 | |
| | 729 | |
| | 730 | FitMomentumResolution(); |
| | 731 | |
| | 732 | } |
| | 733 | |
| | 734 | void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) |
| | 735 | { |
| | 736 | |
| | 737 | // handle and reset the output jet branch |
| | 738 | |
| | 739 | if(fTCAJetsOut)fTCAJetsOut->Delete(); |
| | 740 | if(fTCAJetsOutRan)fTCAJetsOutRan->Delete(); |
| | 741 | if(fTCARandomConesOut)fTCARandomConesOut->Delete(); |
| | 742 | if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete(); |
| | 743 | if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset(); |
| | 744 | |
| | 745 | AliAODJetEventBackground* externalBackground = 0; |
| | 746 | if(!externalBackground&&fBackgroundBranch.Length()){ |
| | 747 | externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data())); |
| | 748 | if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data())); |
| | 749 | if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());; |
| | 750 | } |
| | 751 | // |
| | 752 | // Execute analysis for current event |
| | 753 | // |
| | 754 | AliESDEvent *fESD = 0; |
| | 755 | if(fUseAODTrackInput){ |
| | 756 | fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); |
| | 757 | if(!fAOD){ |
| | 758 | Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput); |
| | 759 | return; |
| | 760 | } |
| | 761 | // fetch the header |
| | 762 | } |
| | 763 | else{ |
| | 764 | // assume that the AOD is in the general output... |
| | 765 | fAOD = AODEvent(); |
| | 766 | if(!fAOD){ |
| | 767 | Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__); |
| | 768 | return; |
| | 769 | } |
| | 770 | if(fDebug>0){ |
| | 771 | fESD = dynamic_cast<AliESDEvent*> (InputEvent()); |
| | 772 | } |
| | 773 | } |
| | 774 | |
| | 775 | //Check if information is provided detector level effects |
| | 776 | if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE; |
| | 777 | if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = kFALSE; |
| | 778 | |
| | 779 | Bool_t selectEvent = false; |
| | 780 | Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB; |
| | 781 | |
| | 782 | Float_t cent = 0; |
| | 783 | Float_t zVtx = 0; |
| | 784 | Int_t cenClass = -1; |
| | 785 | if(fAOD){ |
| | 786 | const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex(); |
| | 787 | TString vtxTitle(vtxAOD->GetTitle()); |
| | 788 | zVtx = vtxAOD->GetZ(); |
| | 789 | |
| | 790 | cent = fAOD->GetHeader()->GetCentrality(); |
| | 791 | if(cent<10)cenClass = 0; |
| | 792 | else if(cent<30)cenClass = 1; |
| | 793 | else if(cent<50)cenClass = 2; |
| | 794 | else if(cent<80)cenClass = 3; |
| | 795 | if(physicsSelection){ |
| | 796 | fh1CentralityPhySel->Fill(cent); |
| | 797 | fh1ZPhySel->Fill(zVtx); |
| | 798 | } |
| | 799 | |
| | 800 | if(fEventSelection){ |
| | 801 | if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){ |
| | 802 | Float_t yvtx = vtxAOD->GetY(); |
| | 803 | Float_t xvtx = vtxAOD->GetX(); |
| | 804 | Float_t r2 = yvtx*yvtx+xvtx*xvtx; |
| | 805 | if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on |
| | 806 | if(physicsSelection){ |
| | 807 | selectEvent = true; |
| | 808 | } |
| | 809 | } |
| | 810 | } |
| | 811 | if(fCentCutUp>0){ |
| | 812 | if(cent<fCentCutLo||cent>fCentCutUp){ |
| | 813 | selectEvent = false; |
| | 814 | } |
| | 815 | } |
| | 816 | }else{ |
| | 817 | selectEvent = true; |
| | 818 | } |
| | 819 | } |
| | 820 | |
| | 821 | |
| | 822 | if(!selectEvent){ |
| | 823 | PostData(1, fHistList); |
| | 824 | return; |
| | 825 | } |
| | 826 | fh1Centrality->Fill(cent); |
| | 827 | fh1Z->Fill(zVtx); |
| | 828 | fh1Trials->Fill("#sum{ntrials}",1); |
| | 829 | |
| | 830 | |
| | 831 | if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__); |
| | 832 | |
| | 833 | // ==== General variables needed |
| | 834 | |
| | 835 | |
| | 836 | |
| | 837 | // we simply fetch the tracks/mc particles as a list of AliVParticles |
| | 838 | |
| | 839 | TList recParticles; |
| | 840 | TList genParticles; |
| | 841 | |
| | 842 | Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec); |
| | 843 | Float_t nCh = recParticles.GetEntries(); |
| | 844 | fh1Nch->Fill(nCh); |
| | 845 | if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries()); |
| | 846 | nT = GetListOfTracks(&genParticles,fTrackTypeGen); |
| | 847 | if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries()); |
| | 848 | |
| | 849 | // find the jets.... |
| | 850 | |
| | 851 | vector<fastjet::PseudoJet> inputParticlesRec; |
| | 852 | vector<fastjet::PseudoJet> inputParticlesRecRan; |
| | 853 | |
| | 854 | // Generate the random cones |
| | 855 | |
| | 856 | AliAODJet vTmpRan(1,0,0,1); |
| | 857 | for(int i = 0; i < recParticles.GetEntries(); i++){ |
| | 858 | AliVParticle *vp = (AliVParticle*)recParticles.At(i); |
| | 859 | |
| | 860 | // Carefull energy is not well determined in real data, should not matter for p_T scheme? |
| | 861 | // we take total momentum here |
| | 862 | |
| | 863 | if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) { |
| | 864 | //Add particles to fastjet in case we are not running toy model |
| | 865 | fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P()); |
| | 866 | jInp.set_user_index(i); |
| | 867 | inputParticlesRec.push_back(jInp); |
| | 868 | } |
| | 869 | else if(fUseDiceEfficiency) { |
| | 870 | |
| | 871 | // Dice to decide if particle is kept or not - toy model for efficiency |
| | 872 | // |
| | 873 | Double_t rnd = fRandom->Uniform(1.); |
| | 874 | Double_t pT = vp->Pt(); |
| | 875 | Double_t eff[3] = {0.}; |
| | 876 | Double_t pTtmp = pT; |
| | 877 | if(pT>10.) pTtmp = 10.; |
| | 878 | Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp)); |
| | 879 | Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp)); |
| | 880 | Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp)); |
| | 881 | Int_t cat[3] = {0}; |
| | 882 | //Sort efficiencies from large to small |
| | 883 | if(eff1>eff2 && eff1>eff3) { |
| | 884 | eff[0] = eff1; |
| | 885 | cat[0] = 1; |
| | 886 | if(eff2>eff3) { |
| | 887 | eff[1] = eff2; |
| | 888 | eff[2] = eff3; |
| | 889 | cat[1] = 2; |
| | 890 | cat[2] = 3; |
| | 891 | } else { |
| | 892 | eff[1] = eff3; |
| | 893 | eff[2] = eff2; |
| | 894 | cat[1] = 3; |
| | 895 | cat[2] = 2; |
| | 896 | } |
| | 897 | } |
| | 898 | else if(eff2>eff1 && eff2>eff3) { |
| | 899 | eff[0] = eff2; |
| | 900 | cat[0] = 2; |
| | 901 | if(eff1>eff3) { |
| | 902 | eff[1] = eff1; |
| | 903 | eff[2] = eff3; |
| | 904 | cat[1] = 1; |
| | 905 | cat[2] = 3; |
| | 906 | } else { |
| | 907 | eff[1] = eff3; |
| | 908 | eff[2] = eff1; |
| | 909 | cat[1] = 3; |
| | 910 | cat[2] = 1; |
| | 911 | } |
| | 912 | } |
| | 913 | else if(eff3>eff1 && eff3>eff2) { |
| | 914 | eff[0] = eff3; |
| | 915 | cat[0] = 3; |
| | 916 | if(eff1>eff2) { |
| | 917 | eff[1] = eff1; |
| | 918 | eff[2] = eff2; |
| | 919 | cat[1] = 1; |
| | 920 | cat[2] = 2; |
| | 921 | } else { |
| | 922 | eff[1] = eff2; |
| | 923 | eff[2] = eff1; |
| | 924 | cat[1] = 2; |
| | 925 | cat[2] = 1; |
| | 926 | } |
| | 927 | } |
| | 928 | |
| | 929 | Double_t sumEff = eff[0]+eff[1]+eff[2]; |
| | 930 | fp1Efficiency->Fill(vp->Pt(),sumEff); |
| | 931 | if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue; |
| | 932 | |
| | 933 | if(fUseTrPtResolutionSmearing) { |
| | 934 | //Smear momentum of generated particle |
| | 935 | Double_t smear = 1.; |
| | 936 | //Select hybrid track category |
| | 937 | if(rnd<=eff[2]) |
| | 938 | smear = GetMomentumSmearing(cat[2],pT); |
| | 939 | else if(rnd<=(eff[2]+eff[1])) |
| | 940 | smear = GetMomentumSmearing(cat[1],pT); |
| | 941 | else |
| | 942 | smear = GetMomentumSmearing(cat[0],pT); |
| | 943 | |
| | 944 | fp1PtResolution->Fill(vp->Pt(),smear); |
| | 945 | |
| | 946 | Double_t sigma = vp->Pt()*smear; |
| | 947 | Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma); |
| | 948 | |
| | 949 | Double_t phi = vp->Phi(); |
| | 950 | Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta()))); |
| | 951 | Double_t pX = pTrec * TMath::Cos(phi); |
| | 952 | Double_t pY = pTrec * TMath::Sin(phi); |
| | 953 | Double_t pZ = pTrec/TMath::Tan(theta); |
| | 954 | Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ); |
| | 955 | |
| | 956 | fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec); |
| | 957 | |
| | 958 | fastjet::PseudoJet jInp(pX,pY,pZ,p); |
| | 959 | jInp.set_user_index(i); |
| | 960 | inputParticlesRec.push_back(jInp); |
| | 961 | |
| | 962 | } |
| | 963 | else { |
| | 964 | fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P()); |
| | 965 | jInp.set_user_index(i); |
| | 966 | inputParticlesRec.push_back(jInp); |
| | 967 | |
| | 968 | } |
| | 969 | |
| | 970 | } |
| | 971 | |
| | 972 | // the randomized input changes eta and phi, but keeps the p_T |
| | 973 | if(i>=fNSkipLeadingRan){// eventually skip the leading particles |
| | 974 | Double_t pT = vp->Pt(); |
| | 975 | Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow; |
| | 976 | Double_t phi = 2.* TMath::Pi() * fRandom->Rndm(); |
| | 977 | |
| | 978 | Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta)); |
| | 979 | Double_t pZ = pT/TMath::Tan(theta); |
| | 980 | |
| | 981 | Double_t pX = pT * TMath::Cos(phi); |
| | 982 | Double_t pY = pT * TMath::Sin(phi); |
| | 983 | Double_t p = TMath::Sqrt(pT*pT+pZ*pZ); |
| | 984 | fastjet::PseudoJet jInpRan(pX,pY,pZ,p); |
| | 985 | |
| | 986 | jInpRan.set_user_index(i); |
| | 987 | inputParticlesRecRan.push_back(jInpRan); |
| | 988 | vTmpRan.SetPxPyPzE(pX,pY,pZ,p); |
| | 989 | } |
| | 990 | |
| | 991 | // fill the tref array, only needed when we write out jets |
| | 992 | if(fTCAJetsOut){ |
| | 993 | if(i == 0){ |
| | 994 | fRef->Delete(); // make sure to delete before placement new... |
| | 995 | if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) { |
| | 996 | new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ... |
| | 997 | } |
| | 998 | } |
| | 999 | if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ... |
| | 1000 | } |
| | 1001 | }// recparticles |
| | 1002 | |
| | 1003 | if(inputParticlesRec.size()==0){ |
| | 1004 | if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__); |
| | 1005 | PostData(1, fHistList); |
| | 1006 | return; |
| | 1007 | } |
| | 1008 | |
| | 1009 | // run fast jet |
| | 1010 | // employ setters for these... |
| | 1011 | |
| | 1012 | |
| | 1013 | // now create the object that holds info about ghosts |
| | 1014 | /* |
| | 1015 | if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){ |
| | 1016 | // reduce CPU time... |
| | 1017 | fGhostArea = 0.5; |
| | 1018 | fActiveAreaRepeats = 0; |
| | 1019 | } |
| | 1020 | */ |
| | 1021 | |
| | 1022 | fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea); |
| | 1023 | fastjet::AreaType areaType = fastjet::active_area; |
| | 1024 | fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec); |
| | 1025 | fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy); |
| | 1026 | fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef); |
| | 1027 | |
| | 1028 | //range where to compute background |
| | 1029 | Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0; |
| | 1030 | phiMin = 0; |
| | 1031 | phiMax = 2*TMath::Pi(); |
| | 1032 | rapMax = fGhostEtamax - fRparam; |
| | 1033 | rapMin = - fGhostEtamax + fRparam; |
| | 1034 | fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax); |
| | 1035 | |
| | 1036 | |
| | 1037 | const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets(); |
| | 1038 | const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets); |
| | 1039 | |
| | 1040 | |
| | 1041 | fh1NJetsRec->Fill(sortedJets.size()); |
| | 1042 | |
| | 1043 | // loop over all jets an fill information, first on is the leading jet |
| | 1044 | |
| | 1045 | Int_t nRecOver = inclusiveJets.size(); |
| | 1046 | Int_t nRec = inclusiveJets.size(); |
| | 1047 | if(inclusiveJets.size()>0){ |
| | 1048 | AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E()); |
| | 1049 | Double_t area = clustSeq.area(sortedJets[0]); |
| | 1050 | leadingJet.SetEffArea(area,0); |
| | 1051 | Float_t pt = leadingJet.Pt(); |
| | 1052 | Int_t nAodOutJets = 0; |
| | 1053 | Int_t nAodOutTracks = 0; |
| | 1054 | AliAODJet *aodOutJet = 0; |
| | 1055 | |
| | 1056 | Int_t iCount = 0; |
| | 1057 | for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){ |
| | 1058 | Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i); |
| | 1059 | while(pt<ptCut&&iCount<nRec){ |
| | 1060 | nRecOver--; |
| | 1061 | iCount++; |
| | 1062 | if(iCount<nRec){ |
| | 1063 | pt = sortedJets[iCount].perp(); |
| | 1064 | } |
| | 1065 | } |
| | 1066 | if(nRecOver<=0)break; |
| | 1067 | fh2NRecJetsPt->Fill(ptCut,nRecOver); |
| | 1068 | } |
| | 1069 | Float_t phi = leadingJet.Phi(); |
| | 1070 | if(phi<0)phi+=TMath::Pi()*2.; |
| | 1071 | Float_t eta = leadingJet.Eta(); |
| | 1072 | Float_t pTback = 0; |
| | 1073 | if(externalBackground){ |
| | 1074 | // carefull has to be filled in a task before |
| | 1075 | // todo, ReArrange to the botom |
| | 1076 | pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged(); |
| | 1077 | } |
| | 1078 | pt = leadingJet.Pt() - pTback; |
| | 1079 | // correlation of leading jet with tracks |
| | 1080 | TIterator *recIter = recParticles.MakeIterator(); |
| | 1081 | recIter->Reset(); |
| | 1082 | AliVParticle *tmpRecTrack = 0; |
| | 1083 | while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){ |
| | 1084 | Float_t tmpPt = tmpRecTrack->Pt(); |
| | 1085 | // correlation |
| | 1086 | Float_t tmpPhi = tmpRecTrack->Phi(); |
| | 1087 | if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; |
| | 1088 | Float_t dPhi = phi - tmpPhi; |
| | 1089 | if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); |
| | 1090 | if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); |
| | 1091 | fh2TracksLeadingJetPhiPt->Fill(dPhi,pt); |
| | 1092 | fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt); |
| | 1093 | if(cenClass>=0){ |
| | 1094 | fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt); |
| | 1095 | fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt); |
| | 1096 | } |
| | 1097 | |
| | 1098 | } |
| | 1099 | |
| | 1100 | |
| | 1101 | TLorentzVector vecareab; |
| | 1102 | for(int j = 0; j < nRec;j++){ |
| | 1103 | AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E()); |
| | 1104 | aodOutJet = 0; |
| | 1105 | nAodOutTracks = 0; |
| | 1106 | Float_t tmpPt = tmpRec.Pt(); |
| | 1107 | |
| | 1108 | if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted... |
| | 1109 | aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec); |
| | 1110 | aodOutJet->GetRefTracks()->Clear(); |
| | 1111 | Double_t area1 = clustSeq.area(sortedJets[j]); |
| | 1112 | aodOutJet->SetEffArea(area1,0); |
| | 1113 | fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]); |
| | 1114 | vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e()); |
| | 1115 | aodOutJet->SetVectorAreaCharged(&vecareab); |
| | 1116 | } |
| | 1117 | |
| | 1118 | |
| | 1119 | Float_t tmpPtBack = 0; |
| | 1120 | if(externalBackground){ |
| | 1121 | // carefull has to be filled in a task before |
| | 1122 | // todo, ReArrange to the botom |
| | 1123 | tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged(); |
| | 1124 | } |
| | 1125 | tmpPt = tmpPt - tmpPtBack; |
| | 1126 | if(tmpPt<0)tmpPt = 0; // avoid negative weights... |
| | 1127 | |
| | 1128 | fh1PtJetsRecIn->Fill(tmpPt); |
| | 1129 | // Fill Spectra with constituentsemacs |
| | 1130 | const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]); |
| | 1131 | |
| | 1132 | fh1NConstRec->Fill(constituents.size()); |
| | 1133 | fh2PtNch->Fill(nCh,tmpPt); |
| | 1134 | fh2PtNchN->Fill(nCh,tmpPt,constituents.size()); |
| | 1135 | fh2NConstPt->Fill(tmpPt,constituents.size()); |
| | 1136 | // loop over constiutents and fill spectrum |
| | 1137 | |
| | 1138 | AliVParticle *partLead = 0; |
| | 1139 | Float_t ptLead = -1; |
| | 1140 | |
| | 1141 | for(unsigned int ic = 0; ic < constituents.size();ic++){ |
| | 1142 | AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index()); |
| | 1143 | if(!part) continue; |
| | 1144 | fh1PtJetConstRec->Fill(part->Pt()); |
| | 1145 | if(aodOutJet){ |
| | 1146 | if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index())); |
| | 1147 | if(part->Pt()>fMaxTrackPtInJet){ |
| | 1148 | aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered); |
| | 1149 | } |
| | 1150 | } |
| | 1151 | if(part->Pt()>ptLead){ |
| | 1152 | ptLead = part->Pt(); |
| | 1153 | partLead = part; |
| | 1154 | } |
| | 1155 | if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt()); |
| | 1156 | } |
| | 1157 | |
| | 1158 | AliAODTrack *aodT = 0; |
| | 1159 | if(partLead&&aodOutJet){ |
| | 1160 | //set pT of leading constituent of jet |
| | 1161 | aodOutJet->SetPtLeading(partLead->Pt()); |
| | 1162 | aodT = dynamic_cast<AliAODTrack*>(partLead); |
| | 1163 | if(aodT){ |
| | 1164 | if(aodT->TestFilterBit(fFilterMaskBestPt)){ |
| | 1165 | aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest); |
| | 1166 | } |
| | 1167 | } |
| | 1168 | } |
| | 1169 | |
| | 1170 | // correlation |
| | 1171 | Float_t tmpPhi = tmpRec.Phi(); |
| | 1172 | Float_t tmpEta = tmpRec.Eta(); |
| | 1173 | if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; |
| | 1174 | if(j==0){ |
| | 1175 | fh1PtJetsLeadingRecIn->Fill(tmpPt); |
| | 1176 | fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta); |
| | 1177 | fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt); |
| | 1178 | fh1NConstLeadingRec->Fill(constituents.size()); |
| | 1179 | fh2NConstLeadingPt->Fill(tmpPt,constituents.size()); |
| | 1180 | continue; |
| | 1181 | } |
| | 1182 | fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta); |
| | 1183 | fh2JetEtaPt->Fill(tmpEta,tmpPt); |
| | 1184 | Float_t dPhi = phi - tmpPhi; |
| | 1185 | if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); |
| | 1186 | if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); |
| | 1187 | Float_t dEta = eta - tmpRec.Eta(); |
| | 1188 | fh2JetsLeadingPhiEta->Fill(dPhi,dEta); |
| | 1189 | fh2JetsLeadingPhiPt->Fill(dPhi,pt); |
| | 1190 | if(cenClass>=0){ |
| | 1191 | fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt); |
| | 1192 | fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt); |
| | 1193 | } |
| | 1194 | fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt); |
| | 1195 | }// loop over reconstructed jets |
| | 1196 | delete recIter; |
| | 1197 | |
| | 1198 | |
| | 1199 | |
| | 1200 | // Add the random cones... |
| | 1201 | if(fNRandomCones>0&&fTCARandomConesOut){ |
| | 1202 | // create a random jet within the acceptance |
| | 1203 | Double_t etaMax = fTrackEtaWindow - fRparam; |
| | 1204 | Int_t nCone = 0; |
| | 1205 | Int_t nConeRan = 0; |
| | 1206 | Double_t pTC = 1; // small number |
| | 1207 | for(int ir = 0;ir < fNRandomCones;ir++){ |
| | 1208 | Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax |
| | 1209 | Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi |
| | 1210 | // massless jet |
| | 1211 | Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); |
| | 1212 | Double_t pZC = pTC/TMath::Tan(thetaC); |
| | 1213 | Double_t pXC = pTC * TMath::Cos(phiC); |
| | 1214 | Double_t pYC = pTC * TMath::Sin(phiC); |
| | 1215 | Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); |
| | 1216 | AliAODJet tmpRecC (pXC,pYC,pZC, pC); |
| | 1217 | bool skip = false; |
| | 1218 | for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets |
| | 1219 | AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E()); |
| | 1220 | if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){ |
| | 1221 | skip = true; |
| | 1222 | break; |
| | 1223 | } |
| | 1224 | } |
| | 1225 | // test for overlap with previous cones to avoid double counting |
| | 1226 | for(int iic = 0;iic<ir;iic++){ |
| | 1227 | AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic); |
| | 1228 | if(iicone){ |
| | 1229 | if(iicone->DeltaR(&tmpRecC)<2.*fRparam){ |
| | 1230 | skip = true; |
| | 1231 | break; |
| | 1232 | } |
| | 1233 | } |
| | 1234 | } |
| | 1235 | if(skip)continue; |
| | 1236 | tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below |
| | 1237 | tmpRecC.SetPtLeading(-1.); |
| | 1238 | if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC); |
| | 1239 | if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC); |
| | 1240 | }// loop over random cones creation |
| | 1241 | |
| | 1242 | |
| | 1243 | // loop over the reconstructed particles and add up the pT in the random cones |
| | 1244 | // maybe better to loop over randomized particles not in the real jets... |
| | 1245 | // but this by definition brings dow average energy in the whole event |
| | 1246 | AliAODJet vTmpRanR(1,0,0,1); |
| | 1247 | for(int i = 0; i < recParticles.GetEntries(); i++){ |
| | 1248 | AliVParticle *vp = (AliVParticle*)recParticles.At(i); |
| | 1249 | if(fTCARandomConesOut){ |
| | 1250 | for(int ir = 0;ir < fNRandomCones;ir++){ |
| | 1251 | AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir); |
| | 1252 | if(jC&&jC->DeltaR(vp)<fRparam){ |
| | 1253 | if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered); |
| | 1254 | jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0); |
| | 1255 | if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt()); |
| | 1256 | } |
| | 1257 | } |
| | 1258 | }// add up energy in cone |
| | 1259 | |
| | 1260 | // the randomized input changes eta and phi, but keeps the p_T |
| | 1261 | if(i>=fNSkipLeadingRan){// eventually skip the leading particles |
| | 1262 | Double_t pTR = vp->Pt(); |
| | 1263 | Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow; |
| | 1264 | Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm(); |
| | 1265 | |
| | 1266 | Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR)); |
| | 1267 | Double_t pZR = pTR/TMath::Tan(thetaR); |
| | 1268 | |
| | 1269 | Double_t pXR = pTR * TMath::Cos(phiR); |
| | 1270 | Double_t pYR = pTR * TMath::Sin(phiR); |
| | 1271 | Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR); |
| | 1272 | vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR); |
| | 1273 | if(fTCARandomConesOutRan){ |
| | 1274 | for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){ |
| | 1275 | AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir); |
| | 1276 | if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){ |
| | 1277 | if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered); |
| | 1278 | jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0); |
| | 1279 | if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt()); |
| | 1280 | } |
| | 1281 | } |
| | 1282 | } |
| | 1283 | } |
| | 1284 | }// loop over recparticles |
| | 1285 | |
| | 1286 | Float_t jetArea = fRparam*fRparam*TMath::Pi(); |
| | 1287 | if(fTCARandomConesOut){ |
| | 1288 | for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){ |
| | 1289 | // rescale the momentum vectors for the random cones |
| | 1290 | |
| | 1291 | AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir); |
| | 1292 | if(rC){ |
| | 1293 | Double_t etaC = rC->Eta(); |
| | 1294 | Double_t phiC = rC->Phi(); |
| | 1295 | // massless jet, unit vector |
| | 1296 | pTC = rC->ChargedBgEnergy(); |
| | 1297 | if(pTC<=0)pTC = 0.001; // for almost empty events |
| | 1298 | Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); |
| | 1299 | Double_t pZC = pTC/TMath::Tan(thetaC); |
| | 1300 | Double_t pXC = pTC * TMath::Cos(phiC); |
| | 1301 | Double_t pYC = pTC * TMath::Sin(phiC); |
| | 1302 | Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); |
| | 1303 | rC->SetPxPyPzE(pXC,pYC,pZC, pC); |
| | 1304 | rC->SetBgEnergy(0,0); |
| | 1305 | rC->SetEffArea(jetArea,0); |
| | 1306 | } |
| | 1307 | } |
| | 1308 | } |
| | 1309 | if(fTCARandomConesOutRan){ |
| | 1310 | for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){ |
| | 1311 | AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir); |
| | 1312 | // same wit random |
| | 1313 | if(rC){ |
| | 1314 | Double_t etaC = rC->Eta(); |
| | 1315 | Double_t phiC = rC->Phi(); |
| | 1316 | // massless jet, unit vector |
| | 1317 | pTC = rC->ChargedBgEnergy(); |
| | 1318 | if(pTC<=0)pTC = 0.001;// for almost empty events |
| | 1319 | Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); |
| | 1320 | Double_t pZC = pTC/TMath::Tan(thetaC); |
| | 1321 | Double_t pXC = pTC * TMath::Cos(phiC); |
| | 1322 | Double_t pYC = pTC * TMath::Sin(phiC); |
| | 1323 | Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); |
| | 1324 | rC->SetPxPyPzE(pXC,pYC,pZC, pC); |
| | 1325 | rC->SetBgEnergy(0,0); |
| | 1326 | rC->SetEffArea(jetArea,0); |
| | 1327 | } |
| | 1328 | } |
| | 1329 | } |
| | 1330 | }// if(fNRandomCones |
| | 1331 | |
| | 1332 | //background estimates:all bckg jets(0) & wo the 2 hardest(1) |
| | 1333 | |
| | 1334 | |
| | 1335 | |
| | 1336 | |
| | 1337 | |
| | 1338 | if(fAODJetBackgroundOut){ |
| | 1339 | vector<fastjet::PseudoJet> jets2=sortedJets; |
| | 1340 | if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); |
| | 1341 | Double_t bkg1=0; |
| | 1342 | Double_t sigma1=0.; |
| | 1343 | Double_t meanarea1=0.; |
| | 1344 | Double_t bkg2=0; |
| | 1345 | Double_t sigma2=0.; |
| | 1346 | Double_t meanarea2=0.; |
| | 1347 | |
| | 1348 | clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true); |
| | 1349 | fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1); |
| | 1350 | |
| | 1351 | // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone)); |
| | 1352 | // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone)); |
| | 1353 | |
| | 1354 | clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true); |
| | 1355 | fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2); |
| | 1356 | // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone)); |
| | 1357 | // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone)); |
| | 1358 | |
| | 1359 | } |
| | 1360 | } |
| | 1361 | |
| | 1362 | |
| | 1363 | |
| | 1364 | |
| | 1365 | |
| | 1366 | // fill track information |
| | 1367 | Int_t nTrackOver = recParticles.GetSize(); |
| | 1368 | // do the same for tracks and jets |
| | 1369 | |
| | 1370 | if(nTrackOver>0){ |
| | 1371 | TIterator *recIter = recParticles.MakeIterator(); |
| | 1372 | AliVParticle *tmpRec = (AliVParticle*)(recIter->Next()); |
| | 1373 | Float_t pt = tmpRec->Pt(); |
| | 1374 | |
| | 1375 | // Printf("Leading track p_t %3.3E",pt); |
| | 1376 | for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){ |
| | 1377 | Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i); |
| | 1378 | while(pt<ptCut&&tmpRec){ |
| | 1379 | nTrackOver--; |
| | 1380 | tmpRec = (AliVParticle*)(recIter->Next()); |
| | 1381 | if(tmpRec){ |
| | 1382 | pt = tmpRec->Pt(); |
| | 1383 | } |
| | 1384 | } |
| | 1385 | if(nTrackOver<=0)break; |
| | 1386 | fh2NRecTracksPt->Fill(ptCut,nTrackOver); |
| | 1387 | } |
| | 1388 | |
| | 1389 | recIter->Reset(); |
| | 1390 | AliVParticle *leading = (AliVParticle*)recParticles.At(0); |
| | 1391 | Float_t phi = leading->Phi(); |
| | 1392 | if(phi<0)phi+=TMath::Pi()*2.; |
| | 1393 | Float_t eta = leading->Eta(); |
| | 1394 | pt = leading->Pt(); |
| | 1395 | while((tmpRec = (AliVParticle*)(recIter->Next()))){ |
| | 1396 | Float_t tmpPt = tmpRec->Pt(); |
| | 1397 | Float_t tmpEta = tmpRec->Eta(); |
| | 1398 | fh1PtTracksRecIn->Fill(tmpPt); |
| | 1399 | fh2TrackEtaPt->Fill(tmpEta,tmpPt); |
| | 1400 | if(tmpRec==leading){ |
| | 1401 | fh1PtTracksLeadingRecIn->Fill(tmpPt); |
| | 1402 | fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt); |
| | 1403 | continue; |
| | 1404 | } |
| | 1405 | // correlation |
| | 1406 | Float_t tmpPhi = tmpRec->Phi(); |
| | 1407 | |
| | 1408 | if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; |
| | 1409 | Float_t dPhi = phi - tmpPhi; |
| | 1410 | if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); |
| | 1411 | if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); |
| | 1412 | Float_t dEta = eta - tmpRec->Eta(); |
| | 1413 | fh2TracksLeadingPhiEta->Fill(dPhi,dEta); |
| | 1414 | fh2TracksLeadingPhiPt->Fill(dPhi,pt); |
| | 1415 | fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt); |
| | 1416 | } |
| | 1417 | delete recIter; |
| | 1418 | } |
| | 1419 | |
| | 1420 | // find the random jets |
| | 1421 | |
| | 1422 | fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef); |
| | 1423 | |
| | 1424 | // fill the jet information from random track |
| | 1425 | const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets(); |
| | 1426 | const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan); |
| | 1427 | |
| | 1428 | fh1NJetsRecRan->Fill(sortedJetsRan.size()); |
| | 1429 | |
| | 1430 | // loop over all jets an fill information, first on is the leading jet |
| | 1431 | |
| | 1432 | Int_t nRecOverRan = inclusiveJetsRan.size(); |
| | 1433 | Int_t nRecRan = inclusiveJetsRan.size(); |
| | 1434 | |
| | 1435 | if(inclusiveJetsRan.size()>0){ |
| | 1436 | AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E()); |
| | 1437 | Float_t pt = leadingJet.Pt(); |
| | 1438 | |
| | 1439 | Int_t iCount = 0; |
| | 1440 | TLorentzVector vecarearanb; |
| | 1441 | |
| | 1442 | for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){ |
| | 1443 | Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i); |
| | 1444 | while(pt<ptCut&&iCount<nRecRan){ |
| | 1445 | nRecOverRan--; |
| | 1446 | iCount++; |
| | 1447 | if(iCount<nRecRan){ |
| | 1448 | pt = sortedJetsRan[iCount].perp(); |
| | 1449 | } |
| | 1450 | } |
| | 1451 | if(nRecOverRan<=0)break; |
| | 1452 | fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan); |
| | 1453 | } |
| | 1454 | Float_t phi = leadingJet.Phi(); |
| | 1455 | if(phi<0)phi+=TMath::Pi()*2.; |
| | 1456 | pt = leadingJet.Pt(); |
| | 1457 | |
| | 1458 | // correlation of leading jet with random tracks |
| | 1459 | |
| | 1460 | for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++) |
| | 1461 | { |
| | 1462 | Float_t tmpPt = inputParticlesRecRan[ip].perp(); |
| | 1463 | // correlation |
| | 1464 | Float_t tmpPhi = inputParticlesRecRan[ip].phi(); |
| | 1465 | if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; |
| | 1466 | Float_t dPhi = phi - tmpPhi; |
| | 1467 | if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); |
| | 1468 | if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); |
| | 1469 | fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt); |
| | 1470 | fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt); |
| | 1471 | } |
| | 1472 | |
| | 1473 | Int_t nAodOutJetsRan = 0; |
| | 1474 | AliAODJet *aodOutJetRan = 0; |
| | 1475 | for(int j = 0; j < nRecRan;j++){ |
| | 1476 | AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E()); |
| | 1477 | Float_t tmpPt = tmpRec.Pt(); |
| | 1478 | fh1PtJetsRecInRan->Fill(tmpPt); |
| | 1479 | // Fill Spectra with constituents |
| | 1480 | const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]); |
| | 1481 | fh1NConstRecRan->Fill(constituents.size()); |
| | 1482 | fh2NConstPtRan->Fill(tmpPt,constituents.size()); |
| | 1483 | fh2PtNchRan->Fill(nCh,tmpPt); |
| | 1484 | fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size()); |
| | 1485 | |
| | 1486 | |
| | 1487 | if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){ |
| | 1488 | aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec); |
| | 1489 | Double_t arearan=clustSeqRan.area(sortedJetsRan[j]); |
| | 1490 | aodOutJetRan->GetRefTracks()->Clear(); |
| | 1491 | aodOutJetRan->SetEffArea(arearan,0); |
| | 1492 | fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]); |
| | 1493 | vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e()); |
| | 1494 | aodOutJetRan->SetVectorAreaCharged(&vecarearanb); |
| | 1495 | |
| | 1496 | } |
| | 1497 | |
| | 1498 | // correlation |
| | 1499 | Float_t tmpPhi = tmpRec.Phi(); |
| | 1500 | if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; |
| | 1501 | |
| | 1502 | if(j==0){ |
| | 1503 | fh1PtJetsLeadingRecInRan->Fill(tmpPt); |
| | 1504 | fh1NConstLeadingRecRan->Fill(constituents.size()); |
| | 1505 | fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size()); |
| | 1506 | continue; |
| | 1507 | } |
| | 1508 | } |
| | 1509 | |
| | 1510 | |
| | 1511 | if(fAODJetBackgroundOut){ |
| | 1512 | Double_t bkg3=0.; |
| | 1513 | Double_t sigma3=0.; |
| | 1514 | Double_t meanarea3=0.; |
| | 1515 | clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true); |
| | 1516 | fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3); |
| | 1517 | // float areaRandomCone = rRandomCone2 *TMath::Pi(); |
| | 1518 | /* |
| | 1519 | fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone)); |
| | 1520 | fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone)); |
| | 1521 | */ |
| | 1522 | } |
| | 1523 | |
| | 1524 | |
| | 1525 | |
| | 1526 | } |
| | 1527 | |
| | 1528 | |
| | 1529 | // do the event selection if activated |
| | 1530 | if(fJetTriggerPtCut>0){ |
| | 1531 | bool select = false; |
| | 1532 | Float_t minPt = fJetTriggerPtCut; |
| | 1533 | /* |
| | 1534 | // hard coded for now ... |
| | 1535 | // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3 |
| | 1536 | if(cent<10)minPt = 50; |
| | 1537 | else if(cent<30)minPt = 42; |
| | 1538 | else if(cent<50)minPt = 28; |
| | 1539 | else if(cent<80)minPt = 18; |
| | 1540 | */ |
| | 1541 | float rho = 0; |
| | 1542 | if(externalBackground)rho = externalBackground->GetBackground(2); |
| | 1543 | if(fTCAJetsOut){ |
| | 1544 | for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){ |
| | 1545 | AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i); |
| | 1546 | Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged(); |
| | 1547 | if(ptSub>=minPt){ |
| | 1548 | select = true; |
| | 1549 | break; |
| | 1550 | } |
| | 1551 | } |
| | 1552 | } |
| | 1553 | |
| | 1554 | if(select){ |
| | 1555 | static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); |
| | 1556 | fh1CentralitySelect->Fill(cent); |
| | 1557 | fh1ZSelect->Fill(zVtx); |
| | 1558 | aodH->SetFillAOD(kTRUE); |
| | 1559 | } |
| | 1560 | } |
| | 1561 | if (fDebug > 2){ |
| | 1562 | if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast()); |
| | 1563 | if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast()); |
| | 1564 | if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast()); |
| | 1565 | if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast()); |
| | 1566 | } |
| | 1567 | PostData(1, fHistList); |
| | 1568 | } |
| | 1569 | |
| | 1570 | void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/) |
| | 1571 | { |
| | 1572 | // |
| | 1573 | // Terminate analysis |
| | 1574 | // |
| | 1575 | if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n"); |
| | 1576 | |
| | 1577 | if(fMomResH1Fit) delete fMomResH1Fit; |
| | 1578 | if(fMomResH2Fit) delete fMomResH2Fit; |
| | 1579 | if(fMomResH3Fit) delete fMomResH3Fit; |
| | 1580 | |
| | 1581 | } |
| | 1582 | |
| | 1583 | |
| | 1584 | Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ |
| | 1585 | |
| | 1586 | // |
| | 1587 | // get list of tracks/particles for different types |
| | 1588 | // |
| | 1589 | |
| | 1590 | if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type); |
| | 1591 | |
| | 1592 | Int_t iCount = 0; |
| | 1593 | if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){ |
| | 1594 | if(type!=kTrackAODextraonly) { |
| | 1595 | AliAODEvent *aod = 0; |
| | 1596 | if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent()); |
| | 1597 | else aod = AODEvent(); |
| | 1598 | if(!aod){ |
| | 1599 | if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__); |
| | 1600 | return iCount; |
| | 1601 | } |
| | 1602 | |
| | 1603 | for(int it = 0;it < aod->GetNumberOfTracks();++it){ |
| | 1604 | AliAODTrack *tr = aod->GetTrack(it); |
| | 1605 | Bool_t bGood = false; |
| | 1606 | if(fFilterType == 0)bGood = true; |
| | 1607 | else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal(); |
| | 1608 | else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal(); |
| | 1609 | if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){ |
| | 1610 | if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap()); |
| | 1611 | continue; |
| | 1612 | } |
| | 1613 | if(TMath::Abs(tr->Eta())>fTrackEtaWindow){ |
| | 1614 | if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks()); |
| | 1615 | continue; |
| | 1616 | } |
| | 1617 | if(tr->Pt()<fTrackPtCut){ |
| | 1618 | if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks()); |
| | 1619 | continue; |
| | 1620 | } |
| | 1621 | if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks()); |
| | 1622 | list->Add(tr); |
| | 1623 | iCount++; |
| | 1624 | } |
| | 1625 | } |
| | 1626 | if(type==kTrackAODextra || type==kTrackAODextraonly) { |
| | 1627 | AliAODEvent *aod = 0; |
| | 1628 | if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent()); |
| | 1629 | else aod = AODEvent(); |
| | 1630 | |
| | 1631 | if(!aod){ |
| | 1632 | return iCount; |
| | 1633 | } |
| | 1634 | TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks")); |
| | 1635 | if(!aodExtraTracks)return iCount; |
| | 1636 | for(int it =0; it<aodExtraTracks->GetEntries(); it++) { |
| | 1637 | AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]); |
| | 1638 | if (!track) continue; |
| | 1639 | |
| | 1640 | AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track); |
| | 1641 | if(!trackAOD)continue; |
| | 1642 | Bool_t bGood = false; |
| | 1643 | if(fFilterType == 0)bGood = true; |
| | 1644 | else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal(); |
| | 1645 | else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal(); |
| | 1646 | if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue; |
| | 1647 | if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue; |
| | 1648 | if(trackAOD->Pt()<fTrackPtCut) continue; |
| | 1649 | list->Add(trackAOD); |
| | 1650 | iCount++; |
| | 1651 | } |
| | 1652 | } |
| | 1653 | } |
| | 1654 | else if (type == kTrackKineAll||type == kTrackKineCharged){ |
| | 1655 | AliMCEvent* mcEvent = MCEvent(); |
| | 1656 | if(!mcEvent)return iCount; |
| | 1657 | // we want to have alivpartilces so use get track |
| | 1658 | for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){ |
| | 1659 | if(!mcEvent->IsPhysicalPrimary(it))continue; |
| | 1660 | AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it); |
| | 1661 | if(type == kTrackKineAll){ |
| | 1662 | if(part->Pt()<fTrackPtCut)continue; |
| | 1663 | list->Add(part); |
| | 1664 | iCount++; |
| | 1665 | } |
| | 1666 | else if(type == kTrackKineCharged){ |
| | 1667 | if(part->Particle()->GetPDG()->Charge()==0)continue; |
| | 1668 | if(part->Pt()<fTrackPtCut)continue; |
| | 1669 | list->Add(part); |
| | 1670 | iCount++; |
| | 1671 | } |
| | 1672 | } |
| | 1673 | } |
| | 1674 | else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) { |
| | 1675 | AliAODEvent *aod = 0; |
| | 1676 | if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent()); |
| | 1677 | else aod = AODEvent(); |
| | 1678 | if(!aod)return iCount; |
| | 1679 | TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName())); |
| | 1680 | if(!tca)return iCount; |
| | 1681 | for(int it = 0;it < tca->GetEntriesFast();++it){ |
| | 1682 | AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it)); |
| | 1683 | if(!part->IsPhysicalPrimary())continue; |
| | 1684 | if(type == kTrackAODMCAll){ |
| | 1685 | if(part->Pt()<fTrackPtCut)continue; |
| | 1686 | list->Add(part); |
| | 1687 | iCount++; |
| | 1688 | } |
| | 1689 | else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){ |
| | 1690 | if(part->Charge()==0)continue; |
| | 1691 | if(part->Pt()<fTrackPtCut)continue; |
| | 1692 | if(kTrackAODMCCharged){ |
| | 1693 | list->Add(part); |
| | 1694 | } |
| | 1695 | else { |
| | 1696 | if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue; |
| | 1697 | list->Add(part); |
| | 1698 | } |
| | 1699 | iCount++; |
| | 1700 | } |
| | 1701 | } |
| | 1702 | }// AODMCparticle |
| | 1703 | list->Sort(); |
| | 1704 | return iCount; |
| | 1705 | } |
| | 1706 | |
| | 1707 | void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() { |
| | 1708 | |
| | 1709 | TFile *f = TFile::Open(fPathTrPtResolution.Data()); |
| | 1710 | |
| | 1711 | if(!f)return; |
| | 1712 | |
| | 1713 | TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt"); |
| | 1714 | TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS"); |
| | 1715 | TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD"); |
| | 1716 | |
| | 1717 | SetSmearResolution(kTRUE); |
| | 1718 | SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD); |
| | 1719 | |
| | 1720 | |
| | 1721 | } |
| | 1722 | |
| | 1723 | void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() { |
| | 1724 | |
| | 1725 | TFile *f = TFile::Open(fPathTrEfficiency.Data()); |
| | 1726 | if(!f)return; |
| | 1727 | |
| | 1728 | TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt"); |
| | 1729 | TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS"); |
| | 1730 | TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD"); |
| | 1731 | |
| | 1732 | SetDiceEfficiency(kTRUE); |
| | 1733 | |
| | 1734 | if(fChangeEfficiencyFraction>0.) { |
| | 1735 | |
| | 1736 | TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin"); |
| | 1737 | |
| | 1738 | for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) { |
| | 1739 | Double_t content = hEffPosGlobSt->GetBinContent(bin); |
| | 1740 | hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction); |
| | 1741 | } |
| | 1742 | |
| | 1743 | SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD); |
| | 1744 | |
| | 1745 | } |
| | 1746 | else |
| | 1747 | SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD); |
| | 1748 | |
| | 1749 | } |
| | 1750 | |
| | 1751 | void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) { |
| | 1752 | |
| | 1753 | // |
| | 1754 | // set mom res profiles |
| | 1755 | // |
| | 1756 | |
| | 1757 | if(fMomResH1) delete fMomResH1; |
| | 1758 | if(fMomResH2) delete fMomResH2; |
| | 1759 | if(fMomResH3) delete fMomResH3; |
| | 1760 | |
| | 1761 | fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1"); |
| | 1762 | fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2"); |
| | 1763 | fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3"); |
| | 1764 | |
| | 1765 | } |
| | 1766 | |
| | 1767 | void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) { |
| | 1768 | // |
| | 1769 | // set tracking efficiency histos |
| | 1770 | // |
| | 1771 | |
| | 1772 | fhEffH1 = (TH1*)h1->Clone("fhEffH1"); |
| | 1773 | fhEffH2 = (TH1*)h2->Clone("fhEffH2"); |
| | 1774 | fhEffH3 = (TH1*)h3->Clone("fhEffH3"); |
| | 1775 | } |
| | 1776 | |
| | 1777 | Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) { |
| | 1778 | |
| | 1779 | // |
| | 1780 | // Get smearing on generated momentum |
| | 1781 | // |
| | 1782 | |
| | 1783 | //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt); |
| | 1784 | |
| | 1785 | TProfile *fMomRes = 0x0; |
| | 1786 | if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes"); |
| | 1787 | if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes"); |
| | 1788 | if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes"); |
| | 1789 | |
| | 1790 | if(!fMomRes) { |
| | 1791 | return 0.; |
| | 1792 | } |
| | 1793 | |
| | 1794 | |
| | 1795 | Double_t smear = 0.; |
| | 1796 | |
| | 1797 | if(pt>20.) { |
| | 1798 | if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt); |
| | 1799 | if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt); |
| | 1800 | if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt); |
| | 1801 | } |
| | 1802 | else { |
| | 1803 | |
| | 1804 | Int_t bin = fMomRes->FindBin(pt); |
| | 1805 | |
| | 1806 | smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin)); |
| | 1807 | |
| | 1808 | } |
| | 1809 | |
| | 1810 | if(fMomRes) delete fMomRes; |
| | 1811 | |
| | 1812 | return smear; |
| | 1813 | } |
| | 1814 | |
| | 1815 | void AliAnalysisTaskJetCluster::FitMomentumResolution() { |
| | 1816 | // |
| | 1817 | // Fit linear function on momentum resolution at high pT |
| | 1818 | // |
| | 1819 | |
| | 1820 | if(!fMomResH1Fit && fMomResH1) { |
| | 1821 | fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.); |
| | 1822 | fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.); |
| | 1823 | fMomResH1Fit ->SetRange(5.,100.); |
| | 1824 | } |
| | 1825 | |
| | 1826 | if(!fMomResH2Fit && fMomResH2) { |
| | 1827 | fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.); |
| | 1828 | fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.); |
| | 1829 | fMomResH2Fit ->SetRange(5.,100.); |
| | 1830 | } |
| | 1831 | |
| | 1832 | if(!fMomResH3Fit && fMomResH3) { |
| | 1833 | fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.); |
| | 1834 | fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.); |
| | 1835 | fMomResH3Fit ->SetRange(5.,100.); |
| | 1836 | } |
| | 1837 | |
| | 1838 | } |
| | 1839 | |
| | 1840 | /* |
| | 1841 | Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){ |
| | 1842 | for(int i = 0; i < particles.GetEntries(); i++){ |
| | 1843 | AliVParticle *vp = (AliVParticle*)particles.At(i); |
| | 1844 | // Carefull energy is not well determined in real data, should not matter for p_T scheme? |
| | 1845 | fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E()); |
| | 1846 | jInp.set_user_index(i); |
| | 1847 | inputParticles.push_back(jInp); |
| | 1848 | } |
| | 1849 | |
| | 1850 | return 0; |
| | 1851 | |
| | 1852 | } |
| | 1853 | */ |
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental