| | 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 <TList.h> |
| | 36 | #include <TLorentzVector.h> |
| | 37 | #include <TClonesArray.h> |
| | 38 | #include "TDatabasePDG.h" |
| | 39 | |
| | 40 | #include "AliAnalysisTaskJetCluster.h" |
| | 41 | #include "AliAnalysisManager.h" |
| | 42 | #include "AliJetFinder.h" |
| | 43 | #include "AliJetHeader.h" |
| | 44 | #include "AliJetReader.h" |
| | 45 | #include "AliESDEvent.h" |
| | 46 | #include "AliAODEvent.h" |
| | 47 | #include "AliAODHandler.h" |
| | 48 | #include "AliAODExtension.h" |
| | 49 | #include "AliAODTrack.h" |
| | 50 | #include "AliAODJet.h" |
| | 51 | #include "AliAODMCParticle.h" |
| | 52 | #include "AliMCEventHandler.h" |
| | 53 | #include "AliMCEvent.h" |
| | 54 | #include "AliStack.h" |
| | 55 | #include "AliGenPythiaEventHeader.h" |
| | 56 | #include "AliJetKineReaderHeader.h" |
| | 57 | #include "AliGenCocktailEventHeader.h" |
| | 58 | #include "AliInputEventHandler.h" |
| | 59 | #include "AliAODJetEventBackground.h" |
| | 60 | |
| | 61 | #include "fastjet/PseudoJet.hh" |
| | 62 | #include "fastjet/ClusterSequenceArea.hh" |
| | 63 | #include "fastjet/AreaDefinition.hh" |
| | 64 | #include "fastjet/JetDefinition.hh" |
| | 65 | // get info on how fastjet was configured |
| | 66 | #include "fastjet/config.h" |
| | 67 | |
| | 68 | |
| | 69 | ClassImp(AliAnalysisTaskJetCluster) |
| | 70 | |
| | 71 | AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){ |
| | 72 | delete fRef; |
| | 73 | delete fRandom; |
| | 74 | |
| | 75 | if(fTCAJetsOut)fTCAJetsOut->Delete(); |
| | 76 | delete fTCAJetsOut; |
| | 77 | if(fTCAJetsOutRan)fTCAJetsOutRan->Delete(); |
| | 78 | delete fTCAJetsOutRan; |
| | 79 | if(fTCARandomConesOut)fTCARandomConesOut->Delete(); |
| | 80 | delete fTCARandomConesOut; |
| | 81 | if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete(); |
| | 82 | delete fTCARandomConesOutRan; |
| | 83 | if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset(); |
| | 84 | delete fAODJetBackgroundOut; |
| | 85 | } |
| | 86 | |
| | 87 | AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), |
| | 88 | fAOD(0x0), |
| | 89 | fAODExtension(0x0), |
| | 90 | fRef(new TRefArray), |
| | 91 | fUseAODTrackInput(kFALSE), |
| | 92 | fUseAODMCInput(kFALSE), |
| | 93 | fUseGlobalSelection(kFALSE), |
| | 94 | fUseBackgroundCalc(kFALSE), |
| | 95 | fFilterMask(0), |
| | 96 | fTrackTypeRec(kTrackUndef), |
| | 97 | fTrackTypeGen(kTrackUndef), |
| | 98 | fNSkipLeadingRan(0), |
| | 99 | fNSkipLeadingCone(0), |
| | 100 | fNRandomCones(0), |
| | 101 | fAvgTrials(1), |
| | 102 | fExternalWeight(1), |
| | 103 | fTrackEtaWindow(0.9), |
| | 104 | fRecEtaWindow(0.5), |
| | 105 | fTrackPtCut(0.), |
| | 106 | fJetOutputMinPt(0.150), |
| | 107 | fJetTriggerPtCut(0), |
| | 108 | fCentCutUp(0), |
| | 109 | fCentCutLo(0), |
| | 110 | fNonStdBranch(""), |
| | 111 | fBackgroundBranch(""), |
| | 112 | fNonStdFile(""), |
| | 113 | fRparam(1.0), |
| | 114 | fAlgorithm(fastjet::kt_algorithm), |
| | 115 | fStrategy(fastjet::Best), |
| | 116 | fRecombScheme(fastjet::BIpt_scheme), |
| | 117 | fAreaType(fastjet::active_area), |
| | 118 | fGhostArea(0.01), |
| | 119 | fActiveAreaRepeats(1), |
| | 120 | fGhostEtamax(1.5), |
| | 121 | fTCAJetsOut(0x0), |
| | 122 | fTCAJetsOutRan(0x0), |
| | 123 | fTCARandomConesOut(0x0), |
| | 124 | fTCARandomConesOutRan(0x0), |
| | 125 | fAODJetBackgroundOut(0x0), |
| | 126 | fRandom(0), |
| | 127 | fh1Xsec(0x0), |
| | 128 | fh1Trials(0x0), |
| | 129 | fh1PtHard(0x0), |
| | 130 | fh1PtHardNoW(0x0), |
| | 131 | fh1PtHardTrials(0x0), |
| | 132 | fh1NJetsRec(0x0), |
| | 133 | fh1NConstRec(0x0), |
| | 134 | fh1NConstLeadingRec(0x0), |
| | 135 | fh1PtJetsRecIn(0x0), |
| | 136 | fh1PtJetsLeadingRecIn(0x0), |
| | 137 | fh1PtJetConstRec(0x0), |
| | 138 | fh1PtJetConstLeadingRec(0x0), |
| | 139 | fh1PtTracksRecIn(0x0), |
| | 140 | fh1PtTracksLeadingRecIn(0x0), |
| | 141 | fh1NJetsRecRan(0x0), |
| | 142 | fh1NConstRecRan(0x0), |
| | 143 | fh1PtJetsLeadingRecInRan(0x0), |
| | 144 | fh1NConstLeadingRecRan(0x0), |
| | 145 | fh1PtJetsRecInRan(0x0), |
| | 146 | fh1PtTracksGenIn(0x0), |
| | 147 | fh1Nch(0x0), |
| | 148 | fh1CentralityPhySel(0x0), |
| | 149 | fh1Centrality(0x0), |
| | 150 | fh1CentralitySelect(0x0), |
| | 151 | fh1ZPhySel(0x0), |
| | 152 | fh1Z(0x0), |
| | 153 | fh1ZSelect(0x0), |
| | 154 | fh2NRecJetsPt(0x0), |
| | 155 | fh2NRecTracksPt(0x0), |
| | 156 | fh2NConstPt(0x0), |
| | 157 | fh2NConstLeadingPt(0x0), |
| | 158 | fh2JetPhiEta(0x0), |
| | 159 | fh2LeadingJetPhiEta(0x0), |
| | 160 | fh2JetEtaPt(0x0), |
| | 161 | fh2LeadingJetEtaPt(0x0), |
| | 162 | fh2TrackEtaPt(0x0), |
| | 163 | fh2LeadingTrackEtaPt(0x0), |
| | 164 | fh2JetsLeadingPhiEta(0x0), |
| | 165 | fh2JetsLeadingPhiPt(0x0), |
| | 166 | fh2TracksLeadingPhiEta(0x0), |
| | 167 | fh2TracksLeadingPhiPt(0x0), |
| | 168 | fh2TracksLeadingJetPhiPt(0x0), |
| | 169 | fh2JetsLeadingPhiPtW(0x0), |
| | 170 | fh2TracksLeadingPhiPtW(0x0), |
| | 171 | fh2TracksLeadingJetPhiPtW(0x0), |
| | 172 | fh2NRecJetsPtRan(0x0), |
| | 173 | fh2NConstPtRan(0x0), |
| | 174 | fh2NConstLeadingPtRan(0x0), |
| | 175 | fh2PtNch(0x0), |
| | 176 | fh2PtNchRan(0x0), |
| | 177 | fh2PtNchN(0x0), |
| | 178 | fh2PtNchNRan(0x0), |
| | 179 | fh2TracksLeadingJetPhiPtRan(0x0), |
| | 180 | fh2TracksLeadingJetPhiPtWRan(0x0), |
| | 181 | fHistList(0x0) |
| | 182 | { |
| | 183 | for(int i = 0;i<3;i++){ |
| | 184 | fh1BiARandomCones[i] = 0; |
| | 185 | fh1BiARandomConesRan[i] = 0; |
| | 186 | } |
| | 187 | for(int i = 0;i<kMaxCent;i++){ |
| | 188 | fh2JetsLeadingPhiPtC[i] = 0; |
| | 189 | fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet |
| | 190 | fh2TracksLeadingJetPhiPtC[i] = 0; |
| | 191 | fh2TracksLeadingJetPhiPtWC[i] = 0; |
| | 192 | } |
| | 193 | } |
| | 194 | |
| | 195 | AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): |
| | 196 | AliAnalysisTaskSE(name), |
| | 197 | fAOD(0x0), |
| | 198 | fAODExtension(0x0), |
| | 199 | fRef(new TRefArray), |
| | 200 | fUseAODTrackInput(kFALSE), |
| | 201 | fUseAODMCInput(kFALSE), |
| | 202 | fUseGlobalSelection(kFALSE), |
| | 203 | fUseBackgroundCalc(kFALSE), |
| | 204 | fFilterMask(0), |
| | 205 | fTrackTypeRec(kTrackUndef), |
| | 206 | fTrackTypeGen(kTrackUndef), |
| | 207 | fNSkipLeadingRan(0), |
| | 208 | fNSkipLeadingCone(0), |
| | 209 | fNRandomCones(0), |
| | 210 | fAvgTrials(1), |
| | 211 | fExternalWeight(1), |
| | 212 | fTrackEtaWindow(0.9), |
| | 213 | fRecEtaWindow(0.5), |
| | 214 | fTrackPtCut(0.), |
| | 215 | fJetOutputMinPt(0.150), |
| | 216 | fJetTriggerPtCut(0), |
| | 217 | fCentCutUp(0), |
| | 218 | fCentCutLo(0), |
| | 219 | fNonStdBranch(""), |
| | 220 | fBackgroundBranch(""), |
| | 221 | fNonStdFile(""), |
| | 222 | fRparam(1.0), |
| | 223 | fAlgorithm(fastjet::kt_algorithm), |
| | 224 | fStrategy(fastjet::Best), |
| | 225 | fRecombScheme(fastjet::BIpt_scheme), |
| | 226 | fAreaType(fastjet::active_area), |
| | 227 | fGhostArea(0.01), |
| | 228 | fActiveAreaRepeats(1), |
| | 229 | fGhostEtamax(1.5), |
| | 230 | fTCAJetsOut(0x0), |
| | 231 | fTCAJetsOutRan(0x0), |
| | 232 | fTCARandomConesOut(0x0), |
| | 233 | fTCARandomConesOutRan(0x0), |
| | 234 | fAODJetBackgroundOut(0x0), |
| | 235 | fRandom(0), |
| | 236 | fh1Xsec(0x0), |
| | 237 | fh1Trials(0x0), |
| | 238 | fh1PtHard(0x0), |
| | 239 | fh1PtHardNoW(0x0), |
| | 240 | fh1PtHardTrials(0x0), |
| | 241 | fh1NJetsRec(0x0), |
| | 242 | fh1NConstRec(0x0), |
| | 243 | fh1NConstLeadingRec(0x0), |
| | 244 | fh1PtJetsRecIn(0x0), |
| | 245 | fh1PtJetsLeadingRecIn(0x0), |
| | 246 | fh1PtJetConstRec(0x0), |
| | 247 | fh1PtJetConstLeadingRec(0x0), |
| | 248 | fh1PtTracksRecIn(0x0), |
| | 249 | fh1PtTracksLeadingRecIn(0x0), |
| | 250 | fh1NJetsRecRan(0x0), |
| | 251 | fh1NConstRecRan(0x0), |
| | 252 | fh1PtJetsLeadingRecInRan(0x0), |
| | 253 | fh1NConstLeadingRecRan(0x0), |
| | 254 | fh1PtJetsRecInRan(0x0), |
| | 255 | fh1PtTracksGenIn(0x0), |
| | 256 | fh1Nch(0x0), |
| | 257 | fh1CentralityPhySel(0x0), |
| | 258 | fh1Centrality(0x0), |
| | 259 | fh1CentralitySelect(0x0), |
| | 260 | fh1ZPhySel(0x0), |
| | 261 | fh1Z(0x0), |
| | 262 | fh1ZSelect(0x0), |
| | 263 | fh2NRecJetsPt(0x0), |
| | 264 | fh2NRecTracksPt(0x0), |
| | 265 | fh2NConstPt(0x0), |
| | 266 | fh2NConstLeadingPt(0x0), |
| | 267 | fh2JetPhiEta(0x0), |
| | 268 | fh2LeadingJetPhiEta(0x0), |
| | 269 | fh2JetEtaPt(0x0), |
| | 270 | fh2LeadingJetEtaPt(0x0), |
| | 271 | fh2TrackEtaPt(0x0), |
| | 272 | fh2LeadingTrackEtaPt(0x0), |
| | 273 | fh2JetsLeadingPhiEta(0x0), |
| | 274 | fh2JetsLeadingPhiPt(0x0), |
| | 275 | fh2TracksLeadingPhiEta(0x0), |
| | 276 | fh2TracksLeadingPhiPt(0x0), |
| | 277 | fh2TracksLeadingJetPhiPt(0x0), |
| | 278 | fh2JetsLeadingPhiPtW(0x0), |
| | 279 | fh2TracksLeadingPhiPtW(0x0), |
| | 280 | fh2TracksLeadingJetPhiPtW(0x0), |
| | 281 | fh2NRecJetsPtRan(0x0), |
| | 282 | fh2NConstPtRan(0x0), |
| | 283 | fh2NConstLeadingPtRan(0x0), |
| | 284 | fh2PtNch(0x0), |
| | 285 | fh2PtNchRan(0x0), |
| | 286 | fh2PtNchN(0x0), |
| | 287 | fh2PtNchNRan(0x0), |
| | 288 | fh2TracksLeadingJetPhiPtRan(0x0), |
| | 289 | fh2TracksLeadingJetPhiPtWRan(0x0), |
| | 290 | fHistList(0x0) |
| | 291 | { |
| | 292 | for(int i = 0;i<3;i++){ |
| | 293 | fh1BiARandomCones[i] = 0; |
| | 294 | fh1BiARandomConesRan[i] = 0; |
| | 295 | } |
| | 296 | for(int i = 0;i<kMaxCent;i++){ |
| | 297 | fh2JetsLeadingPhiPtC[i] = 0; |
| | 298 | fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet |
| | 299 | fh2TracksLeadingJetPhiPtC[i] = 0; |
| | 300 | fh2TracksLeadingJetPhiPtWC[i] = 0; |
| | 301 | } |
| | 302 | DefineOutput(1, TList::Class()); |
| | 303 | } |
| | 304 | |
| | 305 | |
| | 306 | |
| | 307 | Bool_t AliAnalysisTaskJetCluster::Notify() |
| | 308 | { |
| | 309 | // |
| | 310 | // Implemented Notify() to read the cross sections |
| | 311 | // and number of trials from pyxsec.root |
| | 312 | // |
| | 313 | return kTRUE; |
| | 314 | } |
| | 315 | |
| | 316 | void AliAnalysisTaskJetCluster::UserCreateOutputObjects() |
| | 317 | { |
| | 318 | |
| | 319 | // |
| | 320 | // Create the output container |
| | 321 | // |
| | 322 | |
| | 323 | fRandom = new TRandom3(0); |
| | 324 | |
| | 325 | |
| | 326 | // Connect the AOD |
| | 327 | |
| | 328 | |
| | 329 | if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n"); |
| | 330 | |
| | 331 | |
| | 332 | |
| | 333 | if(fNonStdBranch.Length()!=0) |
| | 334 | { |
| | 335 | // only create the output branch if we have a name |
| | 336 | // Create a new branch for jets... |
| | 337 | // -> cleared in the UserExec.... |
| | 338 | // here we can also have the case that the brnaches are written to a separate file |
| | 339 | |
| | 340 | fTCAJetsOut = new TClonesArray("AliAODJet", 0); |
| | 341 | fTCAJetsOut->SetName(fNonStdBranch.Data()); |
| | 342 | AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data()); |
| | 343 | |
| | 344 | |
| | 345 | |
| | 346 | fTCAJetsOutRan = new TClonesArray("AliAODJet", 0); |
| | 347 | fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random")); |
| | 348 | AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data()); |
| | 349 | |
| | 350 | if(fUseBackgroundCalc){ |
| | 351 | if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){ |
| | 352 | fAODJetBackgroundOut = new AliAODJetEventBackground(); |
| | 353 | fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())); |
| | 354 | AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data()); |
| | 355 | } |
| | 356 | } |
| | 357 | // create the branch for the random cones with the same R |
| | 358 | TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone); |
| | 359 | |
| | 360 | if(fNRandomCones>0){ |
| | 361 | if(!AODEvent()->FindListObject(cName.Data())){ |
| | 362 | fTCARandomConesOut = new TClonesArray("AliAODJet", 0); |
| | 363 | fTCARandomConesOut->SetName(cName.Data()); |
| | 364 | AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data()); |
| | 365 | } |
| | 366 | // create the branch with the random for the random cones on the random event |
| | 367 | cName = Form("%sRandomCone_random",fNonStdBranch.Data()); |
| | 368 | if(!AODEvent()->FindListObject(cName.Data())){ |
| | 369 | fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0); |
| | 370 | fTCARandomConesOutRan->SetName(cName.Data()); |
| | 371 | AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data()); |
| | 372 | } |
| | 373 | } |
| | 374 | |
| | 375 | if(fNonStdFile.Length()!=0){ |
| | 376 | // |
| | 377 | // case that we have an AOD extension we need to fetch the jets from the extended output |
| | 378 | // we identify the extension aod event by looking for the branchname |
| | 379 | AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); |
| | 380 | // case that we have an AOD extension we need can fetch the background maybe from the extended output |
| | 381 | fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0); |
| | 382 | } |
| | 383 | } |
| | 384 | |
| | 385 | |
| | 386 | if(!fHistList)fHistList = new TList(); |
| | 387 | fHistList->SetOwner(); |
| | 388 | PostData(1, fHistList); // post data in any case once |
| | 389 | |
| | 390 | Bool_t oldStatus = TH1::AddDirectoryStatus(); |
| | 391 | TH1::AddDirectory(kFALSE); |
| | 392 | |
| | 393 | // |
| | 394 | // Histogram |
| | 395 | |
| | 396 | const Int_t nBinPt = 200; |
| | 397 | Double_t binLimitsPt[nBinPt+1]; |
| | 398 | for(Int_t iPt = 0;iPt <= nBinPt;iPt++){ |
| | 399 | if(iPt == 0){ |
| | 400 | binLimitsPt[iPt] = 0.0; |
| | 401 | } |
| | 402 | else {// 1.0 |
| | 403 | binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0; |
| | 404 | } |
| | 405 | } |
| | 406 | |
| | 407 | const Int_t nBinPhi = 90; |
| | 408 | Double_t binLimitsPhi[nBinPhi+1]; |
| | 409 | for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){ |
| | 410 | if(iPhi==0){ |
| | 411 | binLimitsPhi[iPhi] = -1.*TMath::Pi(); |
| | 412 | } |
| | 413 | else{ |
| | 414 | binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2; |
| | 415 | } |
| | 416 | } |
| | 417 | |
| | 418 | |
| | 419 | |
| | 420 | const Int_t nBinEta = 40; |
| | 421 | Double_t binLimitsEta[nBinEta+1]; |
| | 422 | for(Int_t iEta = 0;iEta<=nBinEta;iEta++){ |
| | 423 | if(iEta==0){ |
| | 424 | binLimitsEta[iEta] = -2.0; |
| | 425 | } |
| | 426 | else{ |
| | 427 | binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1; |
| | 428 | } |
| | 429 | } |
| | 430 | |
| | 431 | const int nChMax = 4000; |
| | 432 | |
| | 433 | fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1); |
| | 434 | fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>"); |
| | 435 | |
| | 436 | fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1); |
| | 437 | fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}"); |
| | 438 | |
| | 439 | |
| | 440 | fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5); |
| | 441 | fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5); |
| | 442 | |
| | 443 | fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5); |
| | 444 | fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5); |
| | 445 | fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5); |
| | 446 | fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5); |
| | 447 | |
| | 448 | |
| | 449 | fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt); |
| | 450 | fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt); |
| | 451 | fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt); |
| | 452 | |
| | 453 | fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 454 | fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 455 | fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 456 | fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 457 | fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 458 | fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); |
| | 459 | fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt); |
| | 460 | fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt); |
| | 461 | fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt); |
| | 462 | fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5); |
| | 463 | |
| | 464 | fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5); |
| | 465 | fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5); |
| | 466 | fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5); |
| | 467 | |
| | 468 | fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25); |
| | 469 | fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25); |
| | 470 | fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25); |
| | 471 | |
| | 472 | fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 473 | fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 474 | fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 475 | // |
| | 476 | |
| | 477 | |
| | 478 | fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 479 | fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 480 | fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 481 | fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5); |
| | 482 | |
| | 483 | fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt); |
| | 484 | 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); |
| | 485 | 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); |
| | 486 | 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); |
| | 487 | |
| | 488 | |
| | 489 | |
| | 490 | fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta", |
| | 491 | nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta); |
| | 492 | fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta", |
| | 493 | nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta); |
| | 494 | |
| | 495 | fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}", |
| | 496 | nBinEta,binLimitsEta,nBinPt,binLimitsPt); |
| | 497 | fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}", |
| | 498 | nBinEta,binLimitsEta,nBinPt,binLimitsPt); |
| | 499 | |
| | 500 | fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}", |
| | 501 | nBinEta,binLimitsEta,nBinPt,binLimitsPt); |
| | 502 | fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}", |
| | 503 | nBinEta,binLimitsEta,nBinPt,binLimitsPt); |
| | 504 | |
| | 505 | |
| | 506 | |
| | 507 | fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta", |
| | 508 | nBinPhi,binLimitsPhi,nBinEta,binLimitsEta); |
| | 509 | fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 510 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 511 | fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta", |
| | 512 | nBinPhi,binLimitsPhi,nBinEta,binLimitsEta); |
| | 513 | fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 514 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 515 | fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 516 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 517 | fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 518 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 519 | |
| | 520 | fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 521 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 522 | fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)", |
| | 523 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 524 | |
| | 525 | fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 526 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 527 | fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", |
| | 528 | nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); |
| | 529 | |
| | 530 | |
| | 531 | if(fNRandomCones>0&&fUseBackgroundCalc){ |
| | 532 | for(int i = 0;i<3;i++){ |
| | 533 | fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100); |
| | 534 | fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100); |
| | 535 | } |
| | 536 | } |
| | 537 | |
| | 538 | for(int i = 0;i < kMaxCent;i++){ |
| | 539 | fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1)); |
| | 540 | fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1)); |
| | 541 | fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1)); |
| | 542 | fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1)); |
| | 543 | } |
| | 544 | |
| | 545 | const Int_t saveLevel = 3; // large save level more histos |
| | 546 | if(saveLevel>0){ |
| | 547 | fHistList->Add(fh1Xsec); |
| | 548 | fHistList->Add(fh1Trials); |
| | 549 | |
| | 550 | fHistList->Add(fh1NJetsRec); |
| | 551 | fHistList->Add(fh1NConstRec); |
| | 552 | fHistList->Add(fh1NConstLeadingRec); |
| | 553 | fHistList->Add(fh1PtJetsRecIn); |
| | 554 | fHistList->Add(fh1PtJetsLeadingRecIn); |
| | 555 | fHistList->Add(fh1PtTracksRecIn); |
| | 556 | fHistList->Add(fh1PtTracksLeadingRecIn); |
| | 557 | fHistList->Add(fh1PtJetConstRec); |
| | 558 | fHistList->Add(fh1PtJetConstLeadingRec); |
| | 559 | fHistList->Add(fh1NJetsRecRan); |
| | 560 | fHistList->Add(fh1NConstRecRan); |
| | 561 | fHistList->Add(fh1PtJetsLeadingRecInRan); |
| | 562 | fHistList->Add(fh1NConstLeadingRecRan); |
| | 563 | fHistList->Add(fh1PtJetsRecInRan); |
| | 564 | fHistList->Add(fh1Nch); |
| | 565 | fHistList->Add(fh1Centrality); |
| | 566 | fHistList->Add(fh1CentralitySelect); |
| | 567 | fHistList->Add(fh1CentralityPhySel); |
| | 568 | fHistList->Add(fh1Z); |
| | 569 | fHistList->Add(fh1ZSelect); |
| | 570 | fHistList->Add(fh1ZPhySel); |
| | 571 | if(fNRandomCones>0&&fUseBackgroundCalc){ |
| | 572 | for(int i = 0;i<3;i++){ |
| | 573 | fHistList->Add(fh1BiARandomCones[i]); |
| | 574 | fHistList->Add(fh1BiARandomConesRan[i]); |
| | 575 | } |
| | 576 | } |
| | 577 | for(int i = 0;i < kMaxCent;i++){ |
| | 578 | fHistList->Add(fh2JetsLeadingPhiPtC[i]); |
| | 579 | fHistList->Add(fh2JetsLeadingPhiPtWC[i]); |
| | 580 | fHistList->Add(fh2TracksLeadingJetPhiPtC[i]); |
| | 581 | fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]); |
| | 582 | } |
| | 583 | |
| | 584 | fHistList->Add(fh2NRecJetsPt); |
| | 585 | fHistList->Add(fh2NRecTracksPt); |
| | 586 | fHistList->Add(fh2NConstPt); |
| | 587 | fHistList->Add(fh2NConstLeadingPt); |
| | 588 | fHistList->Add(fh2PtNch); |
| | 589 | fHistList->Add(fh2PtNchRan); |
| | 590 | fHistList->Add(fh2PtNchN); |
| | 591 | fHistList->Add(fh2PtNchNRan); |
| | 592 | fHistList->Add(fh2JetPhiEta); |
| | 593 | fHistList->Add(fh2LeadingJetPhiEta); |
| | 594 | fHistList->Add(fh2JetEtaPt); |
| | 595 | fHistList->Add(fh2LeadingJetEtaPt); |
| | 596 | fHistList->Add(fh2TrackEtaPt); |
| | 597 | fHistList->Add(fh2LeadingTrackEtaPt); |
| | 598 | fHistList->Add(fh2JetsLeadingPhiEta ); |
| | 599 | fHistList->Add(fh2JetsLeadingPhiPt); |
| | 600 | fHistList->Add(fh2TracksLeadingPhiEta); |
| | 601 | fHistList->Add(fh2TracksLeadingPhiPt); |
| | 602 | fHistList->Add(fh2TracksLeadingJetPhiPt); |
| | 603 | fHistList->Add(fh2JetsLeadingPhiPtW); |
| | 604 | fHistList->Add(fh2TracksLeadingPhiPtW); |
| | 605 | fHistList->Add(fh2TracksLeadingJetPhiPtW); |
| | 606 | fHistList->Add(fh2NRecJetsPtRan); |
| | 607 | fHistList->Add(fh2NConstPtRan); |
| | 608 | fHistList->Add(fh2NConstLeadingPtRan); |
| | 609 | fHistList->Add(fh2TracksLeadingJetPhiPtRan); |
| | 610 | fHistList->Add(fh2TracksLeadingJetPhiPtWRan); |
| | 611 | } |
| | 612 | |
| | 613 | // =========== Switch on Sumw2 for all histos =========== |
| | 614 | for (Int_t i=0; i<fHistList->GetEntries(); ++i) { |
| | 615 | TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i)); |
| | 616 | if (h1){ |
| | 617 | h1->Sumw2(); |
| | 618 | continue; |
| | 619 | } |
| | 620 | THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i)); |
| | 621 | if(hn)hn->Sumw2(); |
| | 622 | } |
| | 623 | TH1::AddDirectory(oldStatus); |
| | 624 | } |
| | 625 | |
| | 626 | void AliAnalysisTaskJetCluster::Init() |
| | 627 | { |
| | 628 | // |
| | 629 | // Initialization |
| | 630 | // |
| | 631 | |
| | 632 | if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n"); |
| | 633 | |
| | 634 | } |
| | 635 | |
| | 636 | void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) |
| | 637 | { |
| | 638 | |
| | 639 | if(fUseGlobalSelection){ |
| | 640 | // no selection by the service task, we continue |
| | 641 | if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__); |
| | 642 | PostData(1, fHistList); |
| | 643 | return; |
| | 644 | } |
| | 645 | |
| | 646 | |
| | 647 | |
| | 648 | // handle and reset the output jet branch |
| | 649 | |
| | 650 | if(fTCAJetsOut)fTCAJetsOut->Delete(); |
| | 651 | if(fTCAJetsOutRan)fTCAJetsOutRan->Delete(); |
| | 652 | if(fTCARandomConesOut)fTCARandomConesOut->Delete(); |
| | 653 | if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete(); |
| | 654 | if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset(); |
| | 655 | |
| | 656 | AliAODJetEventBackground* externalBackground = 0; |
| | 657 | if(!externalBackground&&fBackgroundBranch.Length()){ |
| | 658 | externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data())); |
| | 659 | if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data())); |
| | 660 | if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());; |
| | 661 | } |
| | 662 | // |
| | 663 | // Execute analysis for current event |
| | 664 | // |
| | 665 | AliESDEvent *fESD = 0; |
| | 666 | if(fUseAODTrackInput){ |
| | 667 | fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); |
| | 668 | if(!fAOD){ |
| | 669 | Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput); |
| | 670 | return; |
| | 671 | } |
| | 672 | // fethc the header |
| | 673 | } |
| | 674 | else{ |
| | 675 | // assume that the AOD is in the general output... |
| | 676 | fAOD = AODEvent(); |
| | 677 | if(!fAOD){ |
| | 678 | Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__); |
| | 679 | return; |
| | 680 | } |
| | 681 | if(fDebug>0){ |
| | 682 | fESD = dynamic_cast<AliESDEvent*> (InputEvent()); |
| | 683 | } |
| | 684 | } |
| | 685 | |
| | 686 | Bool_t selectEvent = false; |
| | 687 | Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB; |
| | 688 | |
| | 689 | Float_t cent = 0; |
| | 690 | Float_t zVtx = 0; |
| | 691 | Int_t cenClass = -1; |
| | 692 | if(fAOD){ |
| | 693 | const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex(); |
| | 694 | TString vtxTitle(vtxAOD->GetTitle()); |
| | 695 | zVtx = vtxAOD->GetZ(); |
| | 696 | |
| | 697 | cent = fAOD->GetHeader()->GetCentrality(); |
| | 698 | if(cent<10)cenClass = 0; |
| | 699 | else if(cent<30)cenClass = 1; |
| | 700 | else if(cent<50)cenClass = 2; |
| | 701 | else if(cent<80)cenClass = 3; |
| | 702 | if(physicsSelection){ |
| | 703 | fh1CentralityPhySel->Fill(cent); |
| | 704 | fh1ZPhySel->Fill(zVtx); |
| | 705 | } |
| | 706 | |
| | 707 | |
| | 708 | if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){ |
| | 709 | Float_t yvtx = vtxAOD->GetY(); |
| | 710 | Float_t xvtx = vtxAOD->GetX(); |
| | 711 | Float_t r2 = yvtx*yvtx+xvtx*xvtx; |
| | 712 | if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on |
| | 713 | if(physicsSelection){ |
| | 714 | selectEvent = true; |
| | 715 | } |
| | 716 | } |
| | 717 | } |
| | 718 | if(fCentCutUp>0){ |
| | 719 | if(cent<fCentCutLo||cent>fCentCutUp){ |
| | 720 | selectEvent = false; |
| | 721 | } |
| | 722 | } |
| | 723 | |
| | 724 | } |
| | 725 | if(!selectEvent){ |
| | 726 | PostData(1, fHistList); |
| | 727 | return; |
| | 728 | } |
| | 729 | fh1Centrality->Fill(cent); |
| | 730 | fh1Z->Fill(zVtx); |
| | 731 | fh1Trials->Fill("#sum{ntrials}",1); |
| | 732 | |
| | 733 | |
| | 734 | if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__); |
| | 735 | |
| | 736 | // ==== General variables needed |
| | 737 | |
| | 738 | |
| | 739 | |
| | 740 | // we simply fetch the tracks/mc particles as a list of AliVParticles |
| | 741 | |
| | 742 | TList recParticles; |
| | 743 | TList genParticles; |
| | 744 | |
| | 745 | Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec); |
| | 746 | Float_t nCh = recParticles.GetEntries(); |
| | 747 | fh1Nch->Fill(nCh); |
| | 748 | if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries()); |
| | 749 | nT = GetListOfTracks(&genParticles,fTrackTypeGen); |
| | 750 | if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries()); |
| | 751 | |
| | 752 | // find the jets.... |
| | 753 | |
| | 754 | vector<fastjet::PseudoJet> inputParticlesRec; |
| | 755 | vector<fastjet::PseudoJet> inputParticlesRecRan; |
| | 756 | |
| | 757 | // Generate the random cones |
| | 758 | |
| | 759 | AliAODJet vTmpRan(1,0,0,1); |
| | 760 | for(int i = 0; i < recParticles.GetEntries(); i++){ |
| | 761 | AliVParticle *vp = (AliVParticle*)recParticles.At(i); |
| | 762 | // Carefull energy is not well determined in real data, should not matter for p_T scheme? |
| | 763 | // we take total momentum here |
| | 764 | fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P()); |
| | 765 | jInp.set_user_index(i); |
| | 766 | inputParticlesRec.push_back(jInp); |
| | 767 | |
| | 768 | // the randomized input changes eta and phi, but keeps the p_T |
| | 769 | if(i>=fNSkipLeadingRan){// eventually skip the leading particles |
| | 770 | Double_t pT = vp->Pt(); |
| | 771 | Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow; |
| | 772 | Double_t phi = 2.* TMath::Pi() * fRandom->Rndm(); |
| | 773 | |
| | 774 | Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta)); |
| | 775 | Double_t pZ = pT/TMath::Tan(theta); |
| | 776 | |
| | 777 | Double_t pX = pT * TMath::Cos(phi); |
| | 778 | Double_t pY = pT * TMath::Sin(phi); |
| | 779 | Double_t p = TMath::Sqrt(pT*pT+pZ*pZ); |
| | 780 | fastjet::PseudoJet jInpRan(pX,pY,pZ,p); |
| | 781 | |
| | 782 | jInpRan.set_user_index(i); |
| | 783 | inputParticlesRecRan.push_back(jInpRan); |
| | 784 | vTmpRan.SetPxPyPzE(pX,pY,pZ,p); |
| | 785 | } |
| | 786 | |
| | 787 | // fill the tref array, only needed when we write out jets |
| | 788 | if(fTCAJetsOut){ |
| | 789 | if(i == 0){ |
| | 790 | fRef->Delete(); // make sure to delete before placement new... |
| | 791 | new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); |
| | 792 | } |
| | 793 | fRef->Add(vp); |
| | 794 | } |
| | 795 | }// recparticles |
| | 796 | |
| | 797 | if(inputParticlesRec.size()==0){ |
| | 798 | if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__); |
| | 799 | PostData(1, fHistList); |
| | 800 | return; |
| | 801 | } |
| | 802 | |
| | 803 | // run fast jet |
| | 804 | // employ setters for these... |
| | 805 | |
| | 806 | |
| | 807 | // now create the object that holds info about ghosts |
| | 808 | if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){ |
| | 809 | // reduce CPU time... |
| | 810 | fGhostArea = 0.5; |
| | 811 | fActiveAreaRepeats = 0; |
| | 812 | } |
| | 813 | |
| | 814 | fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea); |
| | 815 | fastjet::AreaType areaType = fastjet::active_area; |
| | 816 | fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec); |
| | 817 | fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy); |
| | 818 | vector <fastjet::PseudoJet> inclusiveJets, sortedJets; |
| | 819 | fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef); |
| | 820 | |
| | 821 | //range where to compute background |
| | 822 | Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0; |
| | 823 | phiMin = 0; |
| | 824 | phiMax = 2*TMath::Pi(); |
| | 825 | rapMax = fGhostEtamax - fRparam; |
| | 826 | rapMin = - fGhostEtamax + fRparam; |
| | 827 | fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax); |
| | 828 | |
| | 829 | |
| | 830 | |
| | 831 | inclusiveJets = clustSeq.inclusive_jets(); |
| | 832 | sortedJets = sorted_by_pt(inclusiveJets); |
| | 833 | |
| | 834 | fh1NJetsRec->Fill(sortedJets.size()); |
| | 835 | |
| | 836 | // loop over all jets an fill information, first on is the leading jet |
| | 837 | |
| | 838 | Int_t nRecOver = inclusiveJets.size(); |
| | 839 | Int_t nRec = inclusiveJets.size(); |
| | 840 | if(inclusiveJets.size()>0){ |
| | 841 | AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E()); |
| | 842 | Double_t area = clustSeq.area(sortedJets[0]); |
| | 843 | leadingJet.SetEffArea(area,0); |
| | 844 | Float_t pt = leadingJet.Pt(); |
| | 845 | Int_t nAodOutJets = 0; |
| | 846 | Int_t nAodOutTracks = 0; |
| | 847 | AliAODJet *aodOutJet = 0; |
| | 848 | |
| | 849 | Int_t iCount = 0; |
| | 850 | for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){ |
| | 851 | Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i); |
| | 852 | while(pt<ptCut&&iCount<nRec){ |
| | 853 | nRecOver--; |
| | 854 | iCount++; |
| | 855 | if(iCount<nRec){ |
| | 856 | pt = sortedJets[iCount].perp(); |
| | 857 | } |
| | 858 | } |
| | 859 | if(nRecOver<=0)break; |
| | 860 | fh2NRecJetsPt->Fill(ptCut,nRecOver); |
| | 861 | } |
| | 862 | Float_t phi = leadingJet.Phi(); |
| | 863 | if(phi<0)phi+=TMath::Pi()*2.; |
| | 864 | Float_t eta = leadingJet.Eta(); |
| | 865 | Float_t pTback = 0; |
| | 866 | if(externalBackground){ |
| | 867 | // carefull has to be filled in a task before |
| | 868 | // todo, ReArrange to the botom |
| | 869 | pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged(); |
| | 870 | } |
| | 871 | pt = leadingJet.Pt() - pTback; |
| | 872 | // correlation of leading jet with tracks |
| | 873 | TIterator *recIter = recParticles.MakeIterator(); |
| | 874 | recIter->Reset(); |
| | 875 | AliVParticle *tmpRecTrack = 0; |
| | 876 | while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){ |
| | 877 | Float_t tmpPt = tmpRecTrack->Pt(); |
| | 878 | // correlation |
| | 879 | Float_t tmpPhi = tmpRecTrack->Phi(); |
| | 880 | if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; |
| | 881 | Float_t dPhi = phi - tmpPhi; |
| | 882 | if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); |
| | 883 | if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); |
| | 884 | fh2TracksLeadingJetPhiPt->Fill(dPhi,pt); |
| | 885 | fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt); |
| | 886 | if(cenClass>=0){ |
| | 887 | fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt); |
| | 888 | fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt); |
| | 889 | } |
| | 890 | |
| | 891 | } |
| | 892 | |
| | 893 | |
| | 894 | TLorentzVector vecareab; |
| | 895 | for(int j = 0; j < nRec;j++){ |
| | 896 | AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E()); |
| | 897 | aodOutJet = 0; |
| | 898 | nAodOutTracks = 0; |
| | 899 | Float_t tmpPt = tmpRec.Pt(); |
| | 900 | |
| | 901 | if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted... |
| | 902 | aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec); |
| | 903 | aodOutJet->GetRefTracks()->Clear(); |
| | 904 | Double_t area1 = clustSeq.area(sortedJets[j]); |
| | 905 | aodOutJet->SetEffArea(area1,0); |
| | 906 | fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]); |
| | 907 | vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e()); |
| | 908 | aodOutJet->SetVectorAreaCharged(&vecareab); |
| | 909 | } |
| | 910 | |
| | 911 | |
| | 912 | Float_t tmpPtBack = 0; |
| | 913 | if(externalBackground){ |
| | 914 | // carefull has to be filled in a task before |
| | 915 | // todo, ReArrange to the botom |
| | 916 | tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged(); |
| | 917 | } |
| | 918 | tmpPt = tmpPt - tmpPtBack; |
| | 919 | if(tmpPt<0)tmpPt = 0; // avoid negative weights... |
| | 920 | |
| | 921 | fh1PtJetsRecIn->Fill(tmpPt); |
| | 922 | // Fill Spectra with constituents |
| | 923 | vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]); |
| | 924 | |
| | 925 | fh1NConstRec->Fill(constituents.size()); |
| | 926 | fh2PtNch->Fill(nCh,tmpPt); |
| | 927 | fh2PtNchN->Fill(nCh,tmpPt,constituents.size()); |
| | 928 | fh2NConstPt->Fill(tmpPt,constituents.size()); |
| | 929 | // loop over constiutents and fill spectrum |
| | 930 | |
| | 931 | for(unsigned int ic = 0; ic < constituents.size();ic++){ |
| | 932 | AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index()); |
| | 933 | fh1PtJetConstRec->Fill(part->Pt()); |
| | 934 | if(aodOutJet){ |
| | 935 | aodOutJet->AddTrack(fRef->At(constituents[ic].user_index())); |
| | 936 | } |
| | 937 | if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt()); |
| | 938 | } |
| | 939 | |
| | 940 | // correlation |
| | 941 | Float_t tmpPhi = tmpRec.Phi(); |
| | 942 | Float_t tmpEta = tmpRec.Eta(); |
| | 943 | if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; |
| | 944 | if(j==0){ |
| | 945 | fh1PtJetsLeadingRecIn->Fill(tmpPt); |
| | 946 | fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta); |
| | 947 | fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt); |
| | 948 | fh1NConstLeadingRec->Fill(constituents.size()); |
| | 949 | fh2NConstLeadingPt->Fill(tmpPt,constituents.size()); |
| | 950 | continue; |
| | 951 | } |
| | 952 | fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta); |
| | 953 | fh2JetEtaPt->Fill(tmpEta,tmpPt); |
| | 954 | Float_t dPhi = phi - tmpPhi; |
| | 955 | if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); |
| | 956 | if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); |
| | 957 | Float_t dEta = eta - tmpRec.Eta(); |
| | 958 | fh2JetsLeadingPhiEta->Fill(dPhi,dEta); |
| | 959 | fh2JetsLeadingPhiPt->Fill(dPhi,pt); |
| | 960 | if(cenClass>=0){ |
| | 961 | fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt); |
| | 962 | fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt); |
| | 963 | } |
| | 964 | fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt); |
| | 965 | }// loop over reconstructed jets |
| | 966 | delete recIter; |
| | 967 | |
| | 968 | |
| | 969 | |
| | 970 | // Add the random cones... |
| | 971 | if(fNRandomCones>0&&fTCARandomConesOut){ |
| | 972 | // create a random jet within the acceptance |
| | 973 | Double_t etaMax = fTrackEtaWindow - fRparam; |
| | 974 | Int_t nCone = 0; |
| | 975 | Int_t nConeRan = 0; |
| | 976 | Double_t pTC = 1; // small number |
| | 977 | for(int ir = 0;ir < fNRandomCones;ir++){ |
| | 978 | Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax |
| | 979 | Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi |
| | 980 | // massless jet |
| | 981 | Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); |
| | 982 | Double_t pZC = pTC/TMath::Tan(thetaC); |
| | 983 | Double_t pXC = pTC * TMath::Cos(phiC); |
| | 984 | Double_t pYC = pTC * TMath::Sin(phiC); |
| | 985 | Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); |
| | 986 | AliAODJet tmpRecC (pXC,pYC,pZC, pC); |
| | 987 | bool skip = false; |
| | 988 | for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets |
| | 989 | AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E()); |
| | 990 | if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){ |
| | 991 | skip = true; |
| | 992 | break; |
| | 993 | } |
| | 994 | } |
| | 995 | // test for overlap with previous cones to avoid double counting |
| | 996 | for(int iic = 0;iic<ir;iic++){ |
| | 997 | AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic); |
| | 998 | if(iicone){ |
| | 999 | if(iicone->DeltaR(&tmpRecC)<2.*fRparam){ |
| | 1000 | skip = true; |
| | 1001 | break; |
| | 1002 | } |
| | 1003 | } |
| | 1004 | } |
| | 1005 | if(skip)continue; |
| | 1006 | tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below |
| | 1007 | if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC); |
| | 1008 | if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC); |
| | 1009 | }// loop over random cones creation |
| | 1010 | |
| | 1011 | |
| | 1012 | // loop over the reconstructed particles and add up the pT in the random cones |
| | 1013 | // maybe better to loop over randomized particles not in the real jets... |
| | 1014 | // but this by definition brings dow average energy in the whole event |
| | 1015 | AliAODJet vTmpRanR(1,0,0,1); |
| | 1016 | for(int i = 0; i < recParticles.GetEntries(); i++){ |
| | 1017 | AliVParticle *vp = (AliVParticle*)recParticles.At(i); |
| | 1018 | if(fTCARandomConesOut){ |
| | 1019 | for(int ir = 0;ir < fNRandomCones;ir++){ |
| | 1020 | AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir); |
| | 1021 | if(jC&&jC->DeltaR(vp)<fRparam){ |
| | 1022 | jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0); |
| | 1023 | } |
| | 1024 | } |
| | 1025 | }// add up energy in cone |
| | 1026 | |
| | 1027 | // the randomized input changes eta and phi, but keeps the p_T |
| | 1028 | if(i>=fNSkipLeadingRan){// eventually skip the leading particles |
| | 1029 | Double_t pTR = vp->Pt(); |
| | 1030 | Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow; |
| | 1031 | Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm(); |
| | 1032 | |
| | 1033 | Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR)); |
| | 1034 | Double_t pZR = pTR/TMath::Tan(thetaR); |
| | 1035 | |
| | 1036 | Double_t pXR = pTR * TMath::Cos(phiR); |
| | 1037 | Double_t pYR = pTR * TMath::Sin(phiR); |
| | 1038 | Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR); |
| | 1039 | vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR); |
| | 1040 | if(fTCARandomConesOutRan){ |
| | 1041 | for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){ |
| | 1042 | AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir); |
| | 1043 | if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){ |
| | 1044 | jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0); |
| | 1045 | } |
| | 1046 | } |
| | 1047 | } |
| | 1048 | } |
| | 1049 | }// loop over recparticles |
| | 1050 | |
| | 1051 | Float_t jetArea = fRparam*fRparam*TMath::Pi(); |
| | 1052 | if(fTCARandomConesOut){ |
| | 1053 | for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){ |
| | 1054 | // rescale the momntum vectors for the random cones |
| | 1055 | |
| | 1056 | AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir); |
| | 1057 | if(rC){ |
| | 1058 | Double_t etaC = rC->Eta(); |
| | 1059 | Double_t phiC = rC->Phi(); |
| | 1060 | // massless jet, unit vector |
| | 1061 | pTC = rC->ChargedBgEnergy(); |
| | 1062 | if(pTC<=0)pTC = 0.001; // for almost empty events |
| | 1063 | Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); |
| | 1064 | Double_t pZC = pTC/TMath::Tan(thetaC); |
| | 1065 | Double_t pXC = pTC * TMath::Cos(phiC); |
| | 1066 | Double_t pYC = pTC * TMath::Sin(phiC); |
| | 1067 | Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); |
| | 1068 | rC->SetPxPyPzE(pXC,pYC,pZC, pC); |
| | 1069 | rC->SetBgEnergy(0,0); |
| | 1070 | rC->SetEffArea(jetArea,0); |
| | 1071 | } |
| | 1072 | } |
| | 1073 | } |
| | 1074 | if(fTCARandomConesOutRan){ |
| | 1075 | for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){ |
| | 1076 | AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir); |
| | 1077 | // same wit random |
| | 1078 | if(rC){ |
| | 1079 | Double_t etaC = rC->Eta(); |
| | 1080 | Double_t phiC = rC->Phi(); |
| | 1081 | // massless jet, unit vector |
| | 1082 | pTC = rC->ChargedBgEnergy(); |
| | 1083 | if(pTC<=0)pTC = 0.001;// for almost empty events |
| | 1084 | Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); |
| | 1085 | Double_t pZC = pTC/TMath::Tan(thetaC); |
| | 1086 | Double_t pXC = pTC * TMath::Cos(phiC); |
| | 1087 | Double_t pYC = pTC * TMath::Sin(phiC); |
| | 1088 | Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); |
| | 1089 | rC->SetPxPyPzE(pXC,pYC,pZC, pC); |
| | 1090 | rC->SetBgEnergy(0,0); |
| | 1091 | rC->SetEffArea(jetArea,0); |
| | 1092 | } |
| | 1093 | } |
| | 1094 | } |
| | 1095 | }// if(fNRandomCones |
| | 1096 | |
| | 1097 | //background estimates:all bckg jets(0) & wo the 2 hardest(1) |
| | 1098 | |
| | 1099 | |
| | 1100 | |
| | 1101 | |
| | 1102 | |
| | 1103 | if(fAODJetBackgroundOut){ |
| | 1104 | vector<fastjet::PseudoJet> jets2=sortedJets; |
| | 1105 | if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); |
| | 1106 | Double_t bkg1=0; |
| | 1107 | Double_t sigma1=0.; |
| | 1108 | Double_t meanarea1=0.; |
| | 1109 | Double_t bkg2=0; |
| | 1110 | Double_t sigma2=0.; |
| | 1111 | Double_t meanarea2=0.; |
| | 1112 | |
| | 1113 | clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true); |
| | 1114 | fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1); |
| | 1115 | |
| | 1116 | // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone)); |
| | 1117 | // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone)); |
| | 1118 | |
| | 1119 | clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true); |
| | 1120 | fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2); |
| | 1121 | // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone)); |
| | 1122 | // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone)); |
| | 1123 | |
| | 1124 | } |
| | 1125 | } |
| | 1126 | |
| | 1127 | |
| | 1128 | |
| | 1129 | |
| | 1130 | |
| | 1131 | // fill track information |
| | 1132 | Int_t nTrackOver = recParticles.GetSize(); |
| | 1133 | // do the same for tracks and jets |
| | 1134 | |
| | 1135 | if(nTrackOver>0){ |
| | 1136 | TIterator *recIter = recParticles.MakeIterator(); |
| | 1137 | AliVParticle *tmpRec = (AliVParticle*)(recIter->Next()); |
| | 1138 | Float_t pt = tmpRec->Pt(); |
| | 1139 | |
| | 1140 | // Printf("Leading track p_t %3.3E",pt); |
| | 1141 | for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){ |
| | 1142 | Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i); |
| | 1143 | while(pt<ptCut&&tmpRec){ |
| | 1144 | nTrackOver--; |
| | 1145 | tmpRec = (AliVParticle*)(recIter->Next()); |
| | 1146 | if(tmpRec){ |
| | 1147 | pt = tmpRec->Pt(); |
| | 1148 | } |
| | 1149 | } |
| | 1150 | if(nTrackOver<=0)break; |
| | 1151 | fh2NRecTracksPt->Fill(ptCut,nTrackOver); |
| | 1152 | } |
| | 1153 | |
| | 1154 | recIter->Reset(); |
| | 1155 | AliVParticle *leading = (AliVParticle*)recParticles.At(0); |
| | 1156 | Float_t phi = leading->Phi(); |
| | 1157 | if(phi<0)phi+=TMath::Pi()*2.; |
| | 1158 | Float_t eta = leading->Eta(); |
| | 1159 | pt = leading->Pt(); |
| | 1160 | while((tmpRec = (AliVParticle*)(recIter->Next()))){ |
| | 1161 | Float_t tmpPt = tmpRec->Pt(); |
| | 1162 | Float_t tmpEta = tmpRec->Eta(); |
| | 1163 | fh1PtTracksRecIn->Fill(tmpPt); |
| | 1164 | fh2TrackEtaPt->Fill(tmpEta,tmpPt); |
| | 1165 | if(tmpRec==leading){ |
| | 1166 | fh1PtTracksLeadingRecIn->Fill(tmpPt); |
| | 1167 | fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt); |
| | 1168 | continue; |
| | 1169 | } |
| | 1170 | // correlation |
| | 1171 | Float_t tmpPhi = tmpRec->Phi(); |
| | 1172 | |
| | 1173 | if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; |
| | 1174 | Float_t dPhi = phi - tmpPhi; |
| | 1175 | if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); |
| | 1176 | if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); |
| | 1177 | Float_t dEta = eta - tmpRec->Eta(); |
| | 1178 | fh2TracksLeadingPhiEta->Fill(dPhi,dEta); |
| | 1179 | fh2TracksLeadingPhiPt->Fill(dPhi,pt); |
| | 1180 | fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt); |
| | 1181 | } |
| | 1182 | delete recIter; |
| | 1183 | } |
| | 1184 | |
| | 1185 | // find the random jets |
| | 1186 | vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan; |
| | 1187 | fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef); |
| | 1188 | |
| | 1189 | inclusiveJetsRan = clustSeqRan.inclusive_jets(); |
| | 1190 | sortedJetsRan = sorted_by_pt(inclusiveJetsRan); |
| | 1191 | |
| | 1192 | // fill the jet information from random track |
| | 1193 | |
| | 1194 | fh1NJetsRecRan->Fill(sortedJetsRan.size()); |
| | 1195 | |
| | 1196 | // loop over all jets an fill information, first on is the leading jet |
| | 1197 | |
| | 1198 | Int_t nRecOverRan = inclusiveJetsRan.size(); |
| | 1199 | Int_t nRecRan = inclusiveJetsRan.size(); |
| | 1200 | |
| | 1201 | if(inclusiveJetsRan.size()>0){ |
| | 1202 | AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E()); |
| | 1203 | Float_t pt = leadingJet.Pt(); |
| | 1204 | |
| | 1205 | Int_t iCount = 0; |
| | 1206 | TLorentzVector vecarearanb; |
| | 1207 | |
| | 1208 | for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){ |
| | 1209 | Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i); |
| | 1210 | while(pt<ptCut&&iCount<nRecRan){ |
| | 1211 | nRecOverRan--; |
| | 1212 | iCount++; |
| | 1213 | if(iCount<nRecRan){ |
| | 1214 | pt = sortedJetsRan[iCount].perp(); |
| | 1215 | } |
| | 1216 | } |
| | 1217 | if(nRecOverRan<=0)break; |
| | 1218 | fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan); |
| | 1219 | } |
| | 1220 | Float_t phi = leadingJet.Phi(); |
| | 1221 | if(phi<0)phi+=TMath::Pi()*2.; |
| | 1222 | pt = leadingJet.Pt(); |
| | 1223 | |
| | 1224 | // correlation of leading jet with random tracks |
| | 1225 | |
| | 1226 | for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++) |
| | 1227 | { |
| | 1228 | Float_t tmpPt = inputParticlesRecRan[ip].perp(); |
| | 1229 | // correlation |
| | 1230 | Float_t tmpPhi = inputParticlesRecRan[ip].phi(); |
| | 1231 | if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; |
| | 1232 | Float_t dPhi = phi - tmpPhi; |
| | 1233 | if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); |
| | 1234 | if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); |
| | 1235 | fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt); |
| | 1236 | fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt); |
| | 1237 | } |
| | 1238 | |
| | 1239 | Int_t nAodOutJetsRan = 0; |
| | 1240 | AliAODJet *aodOutJetRan = 0; |
| | 1241 | for(int j = 0; j < nRecRan;j++){ |
| | 1242 | AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E()); |
| | 1243 | Float_t tmpPt = tmpRec.Pt(); |
| | 1244 | fh1PtJetsRecInRan->Fill(tmpPt); |
| | 1245 | // Fill Spectra with constituents |
| | 1246 | vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]); |
| | 1247 | fh1NConstRecRan->Fill(constituents.size()); |
| | 1248 | fh2NConstPtRan->Fill(tmpPt,constituents.size()); |
| | 1249 | fh2PtNchRan->Fill(nCh,tmpPt); |
| | 1250 | fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size()); |
| | 1251 | |
| | 1252 | |
| | 1253 | if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){ |
| | 1254 | aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec); |
| | 1255 | Double_t arearan=clustSeqRan.area(sortedJetsRan[j]); |
| | 1256 | aodOutJetRan->GetRefTracks()->Clear(); |
| | 1257 | aodOutJetRan->SetEffArea(arearan,0); |
| | 1258 | fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]); |
| | 1259 | vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e()); |
| | 1260 | aodOutJetRan->SetVectorAreaCharged(&vecarearanb); |
| | 1261 | |
| | 1262 | } |
| | 1263 | |
| | 1264 | // correlation |
| | 1265 | Float_t tmpPhi = tmpRec.Phi(); |
| | 1266 | if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; |
| | 1267 | |
| | 1268 | if(j==0){ |
| | 1269 | fh1PtJetsLeadingRecInRan->Fill(tmpPt); |
| | 1270 | fh1NConstLeadingRecRan->Fill(constituents.size()); |
| | 1271 | fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size()); |
| | 1272 | continue; |
| | 1273 | } |
| | 1274 | } |
| | 1275 | |
| | 1276 | |
| | 1277 | if(fAODJetBackgroundOut){ |
| | 1278 | Double_t bkg3=0.; |
| | 1279 | Double_t sigma3=0.; |
| | 1280 | Double_t meanarea3=0.; |
| | 1281 | clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true); |
| | 1282 | fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3); |
| | 1283 | // float areaRandomCone = rRandomCone2 *TMath::Pi(); |
| | 1284 | /* |
| | 1285 | fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone)); |
| | 1286 | fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone)); |
| | 1287 | */ |
| | 1288 | } |
| | 1289 | |
| | 1290 | |
| | 1291 | |
| | 1292 | } |
| | 1293 | |
| | 1294 | |
| | 1295 | // do the event selection if activated |
| | 1296 | if(fJetTriggerPtCut>0){ |
| | 1297 | bool select = false; |
| | 1298 | Float_t minPt = fJetTriggerPtCut; |
| | 1299 | /* |
| | 1300 | // hard coded for now ... |
| | 1301 | // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3 |
| | 1302 | if(cent<10)minPt = 50; |
| | 1303 | else if(cent<30)minPt = 42; |
| | 1304 | else if(cent<50)minPt = 28; |
| | 1305 | else if(cent<80)minPt = 18; |
| | 1306 | */ |
| | 1307 | float rho = 0; |
| | 1308 | if(externalBackground)rho = externalBackground->GetBackground(2); |
| | 1309 | if(fTCAJetsOut){ |
| | 1310 | for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){ |
| | 1311 | AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i); |
| | 1312 | Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged(); |
| | 1313 | if(ptSub>=minPt){ |
| | 1314 | select = true; |
| | 1315 | break; |
| | 1316 | } |
| | 1317 | } |
| | 1318 | } |
| | 1319 | |
| | 1320 | if(select){ |
| | 1321 | static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); |
| | 1322 | fh1CentralitySelect->Fill(cent); |
| | 1323 | fh1ZSelect->Fill(zVtx); |
| | 1324 | aodH->SetFillAOD(kTRUE); |
| | 1325 | } |
| | 1326 | } |
| | 1327 | if (fDebug > 10)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast()); |
| | 1328 | if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__); |
| | 1329 | PostData(1, fHistList); |
| | 1330 | } |
| | 1331 | |
| | 1332 | void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/) |
| | 1333 | { |
| | 1334 | // Terminate analysis |
| | 1335 | // |
| | 1336 | if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n"); |
| | 1337 | } |
| | 1338 | |
| | 1339 | |
| | 1340 | Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ |
| | 1341 | |
| | 1342 | if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type); |
| | 1343 | |
| | 1344 | Int_t iCount = 0; |
| | 1345 | if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){ |
| | 1346 | if(type!=kTrackAODextraonly) { |
| | 1347 | AliAODEvent *aod = 0; |
| | 1348 | if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent()); |
| | 1349 | else aod = AODEvent(); |
| | 1350 | if(!aod){ |
| | 1351 | if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__); |
| | 1352 | return iCount; |
| | 1353 | } |
| | 1354 | for(int it = 0;it < aod->GetNumberOfTracks();++it){ |
| | 1355 | AliAODTrack *tr = aod->GetTrack(it); |
| | 1356 | if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))){ |
| | 1357 | if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap()); |
| | 1358 | continue; |
| | 1359 | } |
| | 1360 | if(TMath::Abs(tr->Eta())>fTrackEtaWindow){ |
| | 1361 | if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks()); |
| | 1362 | continue; |
| | 1363 | } |
| | 1364 | if(tr->Pt()<fTrackPtCut){ |
| | 1365 | if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks()); |
| | 1366 | continue; |
| | 1367 | } |
| | 1368 | if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks()); |
| | 1369 | list->Add(tr); |
| | 1370 | iCount++; |
| | 1371 | } |
| | 1372 | } |
| | 1373 | if(type==kTrackAODextra || type==kTrackAODextraonly) { |
| | 1374 | AliAODEvent *aod = 0; |
| | 1375 | if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent()); |
| | 1376 | else aod = AODEvent(); |
| | 1377 | |
| | 1378 | if(!aod){ |
| | 1379 | return iCount; |
| | 1380 | } |
| | 1381 | TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks")); |
| | 1382 | if(!aodExtraTracks)return iCount; |
| | 1383 | for(int it =0; it<aodExtraTracks->GetEntries(); it++) { |
| | 1384 | AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]); |
| | 1385 | if (!track) continue; |
| | 1386 | |
| | 1387 | AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track); |
| | 1388 | if(!trackAOD)continue; |
| | 1389 | if(trackAOD->Pt()<fTrackPtCut) continue; |
| | 1390 | list->Add(trackAOD); |
| | 1391 | iCount++; |
| | 1392 | } |
| | 1393 | } |
| | 1394 | } |
| | 1395 | else if (type == kTrackKineAll||type == kTrackKineCharged){ |
| | 1396 | AliMCEvent* mcEvent = MCEvent(); |
| | 1397 | if(!mcEvent)return iCount; |
| | 1398 | // we want to have alivpartilces so use get track |
| | 1399 | for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){ |
| | 1400 | if(!mcEvent->IsPhysicalPrimary(it))continue; |
| | 1401 | AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it); |
| | 1402 | if(type == kTrackKineAll){ |
| | 1403 | if(part->Pt()<fTrackPtCut)continue; |
| | 1404 | list->Add(part); |
| | 1405 | iCount++; |
| | 1406 | } |
| | 1407 | else if(type == kTrackKineCharged){ |
| | 1408 | if(part->Particle()->GetPDG()->Charge()==0)continue; |
| | 1409 | if(part->Pt()<fTrackPtCut)continue; |
| | 1410 | list->Add(part); |
| | 1411 | iCount++; |
| | 1412 | } |
| | 1413 | } |
| | 1414 | } |
| | 1415 | else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) { |
| | 1416 | AliAODEvent *aod = 0; |
| | 1417 | if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent()); |
| | 1418 | else aod = AODEvent(); |
| | 1419 | if(!aod)return iCount; |
| | 1420 | TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName())); |
| | 1421 | if(!tca)return iCount; |
| | 1422 | for(int it = 0;it < tca->GetEntriesFast();++it){ |
| | 1423 | AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it)); |
| | 1424 | if(!part->IsPhysicalPrimary())continue; |
| | 1425 | if(type == kTrackAODMCAll){ |
| | 1426 | if(part->Pt()<fTrackPtCut)continue; |
| | 1427 | list->Add(part); |
| | 1428 | iCount++; |
| | 1429 | } |
| | 1430 | else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){ |
| | 1431 | if(part->Charge()==0)continue; |
| | 1432 | if(part->Pt()<fTrackPtCut)continue; |
| | 1433 | if(kTrackAODMCCharged){ |
| | 1434 | list->Add(part); |
| | 1435 | } |
| | 1436 | else { |
| | 1437 | if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue; |
| | 1438 | list->Add(part); |
| | 1439 | } |
| | 1440 | iCount++; |
| | 1441 | } |
| | 1442 | } |
| | 1443 | }// AODMCparticle |
| | 1444 | list->Sort(); |
| | 1445 | return iCount; |
| | 1446 | } |
| | 1447 | |
| | 1448 | /* |
| | 1449 | Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){ |
| | 1450 | for(int i = 0; i < particles.GetEntries(); i++){ |
| | 1451 | AliVParticle *vp = (AliVParticle*)particles.At(i); |
| | 1452 | // Carefull energy is not well determined in real data, should not matter for p_T scheme? |
| | 1453 | fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E()); |
| | 1454 | jInp.set_user_index(i); |
| | 1455 | inputParticles.push_back(jInp); |
| | 1456 | } |
| | 1457 | |
| | 1458 | return 0; |
| | 1459 | |
| | 1460 | } |
| | 1461 | */ |
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental