]>
Commit | Line | Data |
---|---|---|
f3050824 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id:$ */ | |
6ef5bfa4 | 17 | |
f3050824 | 18 | #include <TROOT.h> |
19 | #include <TSystem.h> | |
20 | #include <TChain.h> | |
21 | #include <TFile.h> | |
22 | #include <TList.h> | |
23 | #include <TH1I.h> | |
24 | #include <TH1F.h> | |
25 | #include <TH2F.h> | |
623ed0a0 | 26 | #include <TProfile.h> |
f3050824 | 27 | #include <TCanvas.h> |
28 | #include <TVector3.h> | |
29 | #include <TLorentzVector.h> | |
30 | #include <TMath.h> | |
6ef5bfa4 | 31 | #include <TTree.h> |
32 | #include <TBranch.h> | |
623ed0a0 | 33 | #include <TRandom.h> |
f3050824 | 34 | |
35 | #include "AliAnalysisTaskUE.h" | |
36 | #include "AliAnalysisManager.h" | |
37 | #include "AliMCEventHandler.h" | |
623ed0a0 | 38 | #include "AliMCEvent.h" |
f3050824 | 39 | #include "AliAODEvent.h" |
40 | #include "AliAODInputHandler.h" | |
41 | #include "AliAODHandler.h" | |
42 | #include "AliStack.h" | |
43 | #include "AliAODJet.h" | |
44 | #include "AliAODTrack.h" | |
623ed0a0 | 45 | #include "AliAODMCParticle.h" |
f3050824 | 46 | |
623ed0a0 | 47 | #include "AliGenPythiaEventHeader.h" |
48 | #include "AliAnalysisHelperJetTasks.h" | |
49 | #include "AliStack.h" | |
f3050824 | 50 | #include "AliLog.h" |
51 | ||
52 | // | |
53 | // Analysis class for Underlying Event studies | |
54 | // | |
55 | // Look for correlations on the tranverse regions to | |
56 | // the leading charged jet | |
57 | // | |
58 | // This class needs as input AOD with track and Jets | |
59 | // the output is a list of histograms | |
60 | // | |
61 | // AOD can be either connected to the InputEventHandler | |
62 | // for a chain of AOD files | |
63 | // or | |
64 | // to the OutputEventHandler | |
65 | // for a chain of ESD files, so this case class should be | |
66 | // in the train after the Jet finder | |
67 | // | |
68 | // Arian.Abrahantes.Quintana@cern.ch | |
69 | // Ernesto.Lopez.Torres@cern.ch | |
70 | // | |
71 | ||
72 | ||
73 | ClassImp( AliAnalysisTaskUE) | |
74 | ||
75 | //////////////////////////////////////////////////////////////////////// | |
76 | ||
77 | ||
78 | //____________________________________________________________________ | |
79 | AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name): | |
6ef5bfa4 | 80 | AliAnalysisTask(name, ""), |
c4396fd4 | 81 | fDebug(kFALSE), |
82 | fDeltaAOD(kFALSE), | |
83 | fDeltaAODBranch(""), | |
84 | fArrayJets(0x0), | |
6ef5bfa4 | 85 | fAOD(0x0), |
86 | fAODjets(0x0), | |
c4396fd4 | 87 | //fArrayJets(0x0), |
88 | //fDeltaAOD(kFALSE), | |
89 | //fDeltaAODBranch("jetsAOD"), | |
6ef5bfa4 | 90 | fListOfHistos(0x0), |
91 | fBinsPtInHist(30), | |
92 | fMinJetPtInHist(0.), | |
93 | fMaxJetPtInHist(300.), | |
623ed0a0 | 94 | fUseMCParticleBranch(kFALSE), |
95 | fConstrainDistance(kTRUE), | |
96 | fMinDistance(0.2), | |
97 | fSimulateChJetPt(kFALSE), | |
98 | fUseAliStack(kTRUE), | |
6ef5bfa4 | 99 | fAnaType(1), |
100 | fRegionType(1), | |
101 | fConeRadius(0.7), | |
623ed0a0 | 102 | fConePosition(1), |
6ef5bfa4 | 103 | fAreaReg(1.5393), // Pi*0.7*0.7 |
104 | fUseChPartJet(kFALSE), | |
105 | fUseSingleCharge(kFALSE), | |
106 | fUsePositiveCharge(kTRUE), | |
107 | fOrdering(1), | |
108 | fFilterBit(0xFF), | |
109 | fJetsOnFly(kFALSE), | |
623ed0a0 | 110 | fChJetPtMin(5.0), |
6ef5bfa4 | 111 | fJet1EtaCut(0.2), |
112 | fJet2DeltaPhiCut(2.616), // 150 degrees | |
113 | fJet2RatioPtCut(0.8), | |
114 | fJet3PtCut(15.), | |
115 | fTrackPtCut(0.), | |
116 | fTrackEtaCut(0.9), | |
c4396fd4 | 117 | fAvgTrials(1), |
6ef5bfa4 | 118 | fhNJets(0x0), |
119 | fhEleadingPt(0x0), | |
120 | fhMinRegPtDist(0x0), | |
121 | fhRegionMultMin(0x0), | |
122 | fhMinRegAvePt(0x0), | |
123 | fhMinRegSumPt(0x0), | |
124 | fhMinRegMaxPtPart(0x0), | |
125 | fhMinRegSumPtvsMult(0x0), | |
de6f8090 | 126 | fhdNdEtaPhiDist(0x0), |
6ef5bfa4 | 127 | fhFullRegPartPtDistVsEt(0x0), |
128 | fhTransRegPartPtDistVsEt(0x0), | |
129 | fhRegionSumPtMaxVsEt(0x0), | |
130 | fhRegionMultMax(0x0), | |
131 | fhRegionMultMaxVsEt(0x0), | |
132 | fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0), | |
133 | fhRegionMultMinVsEt(0x0), | |
134 | fhRegionAveSumPtVsEt(0x0), | |
135 | fhRegionDiffSumPtVsEt(0x0), | |
136 | fhRegionAvePartPtMaxVsEt(0x0), | |
137 | fhRegionAvePartPtMinVsEt(0x0), | |
138 | fhRegionMaxPartPtMaxVsEt(0x0), | |
623ed0a0 | 139 | fh1Xsec(0x0), |
140 | fh1Trials(0x0), | |
6ef5bfa4 | 141 | fSettingsTree(0x0)//, fhValidRegion(0x0) |
f3050824 | 142 | { |
143 | // Default constructor | |
144 | // Define input and output slots here | |
145 | // Input slot #0 works with a TChain | |
146 | DefineInput(0, TChain::Class()); | |
147 | // Output slot #0 writes into a TList container | |
148 | DefineOutput(0, TList::Class()); | |
149 | } | |
150 | ||
623ed0a0 | 151 | //______________________________________________________________ |
152 | Bool_t AliAnalysisTaskUE::Notify() | |
153 | { | |
154 | // | |
155 | // Implemented Notify() to read the cross sections | |
156 | // and number of trials from pyxsec.root | |
519378fb | 157 | // Copy from AliAnalysisTaskJFSystematics |
856f9829 | 158 | fAvgTrials = 1; |
623ed0a0 | 159 | TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree(); |
519378fb | 160 | Float_t xsection = 0; |
161 | Float_t ftrials = 1; | |
623ed0a0 | 162 | if(tree){ |
163 | TFile *curfile = tree->GetCurrentFile(); | |
164 | if (!curfile) { | |
165 | Error("Notify","No current file"); | |
166 | return kFALSE; | |
167 | } | |
168 | if(!fh1Xsec||!fh1Trials){ | |
169 | Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__); | |
170 | return kFALSE; | |
171 | } | |
c4396fd4 | 172 | AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials); |
623ed0a0 | 173 | fh1Xsec->Fill("<#sigma>",xsection); |
c4396fd4 | 174 | // construct average trials |
175 | Float_t nEntries = (Float_t)tree->GetTree()->GetEntries(); | |
176 | if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries; | |
177 | } | |
178 | ||
623ed0a0 | 179 | return kTRUE; |
180 | } | |
181 | ||
f3050824 | 182 | //____________________________________________________________________ |
183 | void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/) | |
184 | { | |
6ef5bfa4 | 185 | // Connect the input data |
f3050824 | 186 | |
6ef5bfa4 | 187 | // We need AOD with tracks and jets. |
188 | // Since AOD can be either connected to the InputEventHandler (input chain fron AOD files) | |
1e996f55 | 189 | // or to the OutputEventHandler ( AOD is create by a previus task in the train ) |
6ef5bfa4 | 190 | // we need to check where it is and get the pointer to AODEvent in the right way |
f3050824 | 191 | |
c4396fd4 | 192 | // Delta AODs are also implemented |
193 | ||
194 | ||
195 | if (fDebug > 1) AliInfo("ConnectInputData() "); | |
f3050824 | 196 | |
197 | TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); | |
198 | ||
c4396fd4 | 199 | if( handler && handler->InheritsFrom("AliAODInputHandler") ) { // input AOD |
200 | fAOD = ((AliAODInputHandler*)handler)->GetEvent(); | |
201 | if (fDebug > 1) AliInfo(" ==== Tracks from AliAODInputHandler"); | |
202 | // Case when jets are reconstructed on the fly from AOD tracks | |
203 | // (the Jet Finder is using the AliJetAODReader) of InputEventHandler | |
204 | // and put in the OutputEventHandler AOD. Useful whe you want to reconstruct jets with | |
205 | // different parameters to default ones stored in the AOD or to use a different algorithm | |
206 | if( fJetsOnFly ) { | |
207 | handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler(); | |
208 | if( handler && handler->InheritsFrom("AliAODHandler") ) { | |
209 | fAODjets = ((AliAODHandler*)handler)->GetAOD(); | |
210 | if (fDebug > 1) AliInfo(" ==== Jets from AliAODHandler (on the fly)"); | |
211 | } | |
212 | } else { | |
213 | fAODjets = fAOD; | |
214 | if (fDebug > 1) AliInfo(" ==== Jets from AliAODInputHandler"); | |
215 | } | |
216 | } else { //output AOD | |
6ef5bfa4 | 217 | handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler(); |
218 | if( handler && handler->InheritsFrom("AliAODHandler") ) { | |
c4396fd4 | 219 | fAOD = ((AliAODHandler*)handler)->GetAOD(); |
6ef5bfa4 | 220 | fAODjets = fAOD; |
c4396fd4 | 221 | if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler"); |
6ef5bfa4 | 222 | } else { |
223 | AliFatal("I can't get any AOD Event Handler"); | |
224 | return; | |
225 | } | |
226 | } | |
f3050824 | 227 | } |
228 | ||
229 | //____________________________________________________________________ | |
230 | void AliAnalysisTaskUE::CreateOutputObjects() | |
231 | { | |
6ef5bfa4 | 232 | // Create the output container |
233 | // | |
f3050824 | 234 | if (fDebug > 1) AliInfo("CreateOutPutData()"); |
235 | // | |
236 | // Histograms | |
1e996f55 | 237 | |
238 | OpenFile(0); | |
f3050824 | 239 | CreateHistos(); |
1e996f55 | 240 | // fListOfHistos->SetOwner(kTRUE); |
241 | ||
242 | ||
f3050824 | 243 | } |
244 | ||
245 | ||
246 | //____________________________________________________________________ | |
247 | void AliAnalysisTaskUE::Exec(Option_t */*option*/) | |
248 | { | |
6ef5bfa4 | 249 | // Execute analysis for current event |
250 | // | |
251 | if ( fDebug > 3 ) AliInfo( " Processing event..." ); | |
519378fb | 252 | // fetch the pythia header info and get the trials |
253 | AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); | |
254 | Float_t nTrials = 1; | |
c4396fd4 | 255 | if (mcHandler) { |
256 | AliMCEvent* mcEvent = mcHandler->MCEvent(); | |
257 | if (mcEvent) { | |
258 | AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent); | |
259 | if(pythiaGenHeader){ | |
260 | nTrials = pythiaGenHeader->Trials(); | |
261 | } | |
519378fb | 262 | } |
263 | } | |
c4396fd4 | 264 | fh1Trials->Fill("#sum{ntrials}",fAvgTrials); |
265 | AnalyseUE(); | |
6ef5bfa4 | 266 | |
f3050824 | 267 | // Post the data |
268 | PostData(0, fListOfHistos); | |
269 | } | |
f3050824 | 270 | |
271 | //____________________________________________________________________ | |
272 | void AliAnalysisTaskUE::AnalyseUE() | |
273 | { | |
274 | // | |
275 | // Look for correlations on the tranverse regions to | |
276 | // the leading charged jet | |
277 | // | |
6ef5bfa4 | 278 | |
f3050824 | 279 | // ------------------------------------------------ |
280 | // Find Leading Jets 1,2,3 | |
281 | // (could be skipped if Jets are sort by Pt...) | |
282 | Double_t maxPtJet1 = 0.; | |
283 | Int_t index1 = -1; | |
284 | Double_t maxPtJet2 = 0.; // jet 2 need for back to back inclusive | |
285 | Int_t index2 = -1; | |
286 | Double_t maxPtJet3 = 0.; // jet 3 need for back to back exclusive | |
287 | Int_t index3 = -1; | |
6ef5bfa4 | 288 | TVector3 jetVect[3]; |
289 | Int_t nJets = 0; | |
290 | ||
c4396fd4 | 291 | |
6ef5bfa4 | 292 | if( !fUseChPartJet ) { |
293 | ||
294 | // Use AOD Jets | |
c4396fd4 | 295 | if(fDeltaAOD){ |
296 | if (fDebug > 1) AliInfo(" ==== Jets From Delta-AODs !"); | |
297 | if (fDebug > 1) AliInfo(Form(" ==== Reading Branch: %s ",fDeltaAODBranch.Data())); | |
298 | fArrayJets = | |
299 | (TClonesArray*)fAODjets->GetList()->FindObject(fDeltaAODBranch.Data()); | |
300 | if (!fArrayJets){ | |
301 | AliFatal(" No jet-array! "); | |
302 | return; | |
303 | } | |
304 | nJets=fArrayJets->GetEntries(); | |
305 | }else{ | |
306 | if (fDebug > 1) AliInfo(" ==== Read Standard-AODs !"); | |
307 | nJets = fAODjets->GetNJets(); | |
308 | } | |
309 | //printf("AOD %d jets \n", nJets); | |
310 | ||
6ef5bfa4 | 311 | for( Int_t i=0; i<nJets; ++i ) { |
c4396fd4 | 312 | AliAODJet* jet; |
313 | if (fDeltaAOD){ | |
314 | jet =(AliAODJet*)fArrayJets->At(i); | |
315 | }else{ | |
316 | jet = fAODjets->GetJet(i); | |
317 | } | |
318 | Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!! | |
6ef5bfa4 | 319 | if( jetPt > maxPtJet1 ) { |
c4396fd4 | 320 | maxPtJet3 = maxPtJet2; index3 = index2; |
321 | maxPtJet2 = maxPtJet1; index2 = index1; | |
322 | maxPtJet1 = jetPt; index1 = i; | |
6ef5bfa4 | 323 | } else if( jetPt > maxPtJet2 ) { |
c4396fd4 | 324 | maxPtJet3 = maxPtJet2; index3 = index2; |
325 | maxPtJet2 = jetPt; index2 = i; | |
6ef5bfa4 | 326 | } else if( jetPt > maxPtJet3 ) { |
c4396fd4 | 327 | maxPtJet3 = jetPt; index3 = i; |
6ef5bfa4 | 328 | } |
329 | } | |
330 | if( index1 != -1 ) { | |
c4396fd4 | 331 | AliAODJet *jet = 0; |
332 | if (fDeltaAOD) { | |
333 | jet =(AliAODJet*) fArrayJets->At(index1); | |
334 | }else{ | |
335 | jet = fAODjets->GetJet(index1); | |
336 | } | |
6ef5bfa4 | 337 | jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); |
338 | } | |
339 | if( index2 != -1 ) { | |
c4396fd4 | 340 | AliAODJet* jet = 0; |
341 | if (fDeltaAOD) { | |
342 | jet= (AliAODJet*) fArrayJets->At(index2); | |
343 | }else{ | |
344 | jet= fAODjets->GetJet(index2); | |
345 | } | |
6ef5bfa4 | 346 | jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); |
347 | } | |
348 | if( index3 != -1 ) { | |
c4396fd4 | 349 | AliAODJet* jet = 0; |
350 | if (fDeltaAOD) { | |
351 | jet= (AliAODJet*) fArrayJets->At(index3); | |
352 | }else{ | |
353 | fAODjets->GetJet(index3); | |
354 | } | |
355 | jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); | |
6ef5bfa4 | 356 | } |
357 | ||
358 | } else { | |
359 | ||
360 | // Use "Charged Particle Jets" | |
361 | TObjArray* jets = FindChargedParticleJets(); | |
362 | if( jets ) { | |
363 | nJets = jets->GetEntriesFast(); | |
364 | if( nJets > 0 ) { | |
c4396fd4 | 365 | index1 = 0; |
366 | AliAODJet* jet = (AliAODJet*)jets->At(0); | |
367 | maxPtJet1 = jet->Pt(); | |
368 | jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); | |
6ef5bfa4 | 369 | } |
370 | if( nJets > 1 ) { | |
c4396fd4 | 371 | index2 = 1; |
372 | AliAODJet* jet = (AliAODJet*)jets->At(1); | |
373 | maxPtJet1 = jet->Pt(); | |
374 | jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); | |
6ef5bfa4 | 375 | } |
376 | if( nJets > 2 ) { | |
c4396fd4 | 377 | index3 = 2; |
378 | AliAODJet* jet = (AliAODJet*)jets->At(2); | |
379 | maxPtJet1 = jet->Pt(); | |
380 | jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); | |
6ef5bfa4 | 381 | } |
382 | ||
383 | jets->Delete(); | |
384 | delete jets; | |
f3050824 | 385 | } |
386 | } | |
6ef5bfa4 | 387 | |
c4396fd4 | 388 | |
389 | fhNJets->Fill(nJets); | |
f3050824 | 390 | |
6ef5bfa4 | 391 | if( fDebug > 1 ) { |
392 | if( index1 < 0 ) { | |
393 | AliInfo("\n Skipping Event, not jet found..."); | |
394 | return; | |
395 | } else | |
396 | AliInfo(Form("\n Pt Leading Jet = %6.1f eta=%5.3f ", maxPtJet1, jetVect[0].Eta() )); | |
397 | } | |
398 | ||
f3050824 | 399 | |
400 | // ---------------------------------------------- | |
6ef5bfa4 | 401 | // Cut events by jets topology |
f3050824 | 402 | // fAnaType: |
6ef5bfa4 | 403 | // 1 = inclusive, |
f3050824 | 404 | // - Jet1 |eta| < fJet1EtaCut |
405 | // 2 = back to back inclusive | |
406 | // - fulfill case 1 | |
407 | // - |Jet1.Phi - Jet2.Phi| > fJet2DeltaPhiCut | |
408 | // - Jet2.Pt/Jet1Pt > fJet2RatioPtCut | |
6ef5bfa4 | 409 | // 3 = back to back exclusive |
f3050824 | 410 | // - fulfill case 2 |
411 | // - Jet3.Pt < fJet3PtCut | |
412 | ||
6ef5bfa4 | 413 | if( index1 < 0 || TMath::Abs(jetVect[0].Eta()) > fJet1EtaCut) { |
414 | if( fDebug > 1 ) AliInfo("\n Skipping Event...Jet1 |eta| > fJet1EtaCut"); | |
f3050824 | 415 | return; |
416 | } | |
f3050824 | 417 | // back to back inclusive |
6ef5bfa4 | 418 | if( fAnaType > 1 && index2 == -1 ) { |
419 | if( fDebug > 1 ) AliInfo("\n Skipping Event... no second Jet found"); | |
420 | return; | |
421 | } | |
422 | if( fAnaType > 1 && index2 > -1 ) { | |
f3050824 | 423 | if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut || |
6ef5bfa4 | 424 | maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) { |
425 | if( fDebug > 1 ) AliInfo("\n Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut"); | |
426 | return; | |
f3050824 | 427 | } |
f3050824 | 428 | } |
f3050824 | 429 | // back to back exclusive |
6ef5bfa4 | 430 | if( fAnaType > 2 && index3 > -1 ) { |
431 | if( maxPtJet3 > fJet3PtCut ) { | |
432 | if( fDebug > 1 ) AliInfo("\n Skipping Event... Jet3.Pt > fJet3PtCut "); | |
f3050824 | 433 | return; |
434 | } | |
f3050824 | 435 | } |
f3050824 | 436 | |
623ed0a0 | 437 | //fhEleadingPt->Fill( maxPtJet1 ); |
6ef5bfa4 | 438 | //Area for Normalization Purpose at Display histos |
439 | SetRegionArea(jetVect); | |
440 | ||
f3050824 | 441 | // ---------------------------------------------- |
442 | // Find max and min regions | |
443 | Double_t sumPtRegionPosit = 0.; | |
444 | Double_t sumPtRegionNegat = 0.; | |
445 | Double_t maxPartPtRegion = 0.; | |
446 | Int_t nTrackRegionPosit = 0; | |
447 | Int_t nTrackRegionNegat = 0; | |
448 | static const Double_t k270rad = 270.*TMath::Pi()/180.; | |
449 | ||
623ed0a0 | 450 | if (!fUseMCParticleBranch){ |
451 | fhEleadingPt->Fill( maxPtJet1 ); | |
452 | Int_t nTracks = fAOD->GetNTracks(); | |
6ef5bfa4 | 453 | |
623ed0a0 | 454 | for (Int_t ipart=0; ipart<nTracks; ++ipart) { |
455 | ||
456 | AliAODTrack* part = fAOD->GetTrack( ipart ); | |
457 | if ( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection | |
458 | if (!part->IsPrimaryCandidate()) continue; // reject whatever is not linked to collision point | |
459 | // PID Selection: Reject everything but hadrons | |
460 | Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion || | |
461 | part->GetMostProbablePID()==AliAODTrack::kKaon || | |
462 | part->GetMostProbablePID()==AliAODTrack::kProton; | |
463 | if ( !isHadron ) continue; | |
464 | ||
465 | if ( !part->Charge() ) continue; //Only charged | |
466 | if ( fUseSingleCharge ) { // Charge selection | |
467 | if ( fUsePositiveCharge && part->Charge() < 0.) continue; // keep Positives | |
468 | if ( !fUsePositiveCharge && part->Charge() > 0.) continue; // keep Negatives | |
469 | } | |
470 | ||
471 | if ( part->Pt() < fTrackPtCut ) continue; | |
472 | if( TMath::Abs(part->Eta()) > fTrackEtaCut ) continue; | |
473 | ||
474 | TVector3 partVect(part->Px(), part->Py(), part->Pz()); | |
475 | ||
476 | Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad; | |
477 | if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi(); | |
478 | fhdNdEtaPhiDist->Fill( deltaPhi ); | |
479 | ||
480 | fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); | |
481 | ||
482 | Int_t region = IsTrackInsideRegion( jetVect, &partVect ); | |
483 | ||
484 | if (region > 0) { | |
485 | if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt(); | |
486 | sumPtRegionPosit += part->Pt(); | |
487 | nTrackRegionPosit++; | |
488 | fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); | |
489 | } | |
490 | if (region < 0) { | |
491 | if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt(); | |
492 | sumPtRegionNegat += part->Pt(); | |
493 | nTrackRegionNegat++; | |
494 | fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); | |
495 | } | |
496 | }//end loop AOD tracks | |
497 | } | |
498 | else { | |
6ef5bfa4 | 499 | |
623ed0a0 | 500 | // this is the part we only use when we have MC information |
501 | // More than a test for values of it also resumes the reconstruction efficiency of jets | |
502 | // As commented bellow if available for the data, we try to pair reconstructed jets with simulated ones | |
503 | // afterwards we kept angular variables of MC jet to perform UE analysis over MC particles | |
504 | // TODO: Handle Multiple jet environment. 06/2009 just suited for inclusive jet condition ( fAnaType = 1 ) | |
6ef5bfa4 | 505 | |
623ed0a0 | 506 | AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); |
507 | if (!mcHandler) { | |
508 | Printf("ERROR: Could not retrieve MC event handler"); | |
509 | return; | |
510 | } | |
511 | ||
512 | AliMCEvent* mcEvent = mcHandler->MCEvent(); | |
513 | if (!mcEvent) { | |
514 | Printf("ERROR: Could not retrieve MC event"); | |
515 | return; | |
516 | } | |
517 | AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent); | |
518 | if(!pythiaGenHeader){ | |
519 | return; | |
520 | } | |
f3050824 | 521 | |
623ed0a0 | 522 | //Get Jets from MC header |
523 | Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets(); | |
524 | AliAODJet pythiaGenJets[4]; | |
525 | TVector3 jetVectnew[4]; | |
526 | Int_t iCount = 0; | |
527 | for(int ip = 0;ip < nPythiaGenJets;++ip){ | |
528 | if (iCount>3) break; | |
529 | Float_t p[4]; | |
530 | pythiaGenHeader->TriggerJet(ip,p); | |
531 | TVector3 tempVect(p[0],p[1],p[2]); | |
532 | if ( TMath::Abs(tempVect.Eta())>fJet1EtaCut ) continue; | |
533 | pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]); | |
534 | jetVectnew[iCount].SetXYZ(pythiaGenJets[iCount].Px(), pythiaGenJets[iCount].Py(), pythiaGenJets[iCount].Pz()); | |
535 | iCount++; | |
536 | } | |
6ef5bfa4 | 537 | |
623ed0a0 | 538 | if (!iCount) return;// no jet in eta acceptance |
f3050824 | 539 | |
623ed0a0 | 540 | //Search the index of the nearest MC jet to the leading jet reconstructed from the input data |
541 | Int_t index = 0; | |
542 | if (fConstrainDistance){ | |
543 | Float_t deltaR = 0.; | |
544 | Float_t dRTemp = 0.; | |
545 | for (Int_t i=0; i<iCount; i++){ | |
546 | if (!i) { | |
547 | dRTemp = jetVectnew[i].DeltaR(jetVect[0]); | |
548 | index = i; | |
549 | } | |
550 | deltaR = jetVectnew[i].DeltaR(jetVect[0]); | |
551 | if (deltaR < dRTemp){ | |
552 | index = i; | |
553 | dRTemp = deltaR; | |
554 | } | |
555 | } | |
556 | ||
557 | if (jetVectnew[index].DeltaR(jetVect[0]) > fMinDistance) return; | |
558 | } | |
559 | //Let's add some taste to jet and simulate pt of charged alone | |
560 | //eta and phi are kept as original | |
561 | //Play a Normal Distribution | |
562 | Float_t random = 1.; | |
563 | if (fSimulateChJetPt){ | |
564 | while(1){ | |
565 | random = gRandom->Gaus(0.6,0.25); | |
566 | if (random > 0. && random < 1. && | |
567 | (random * jetVectnew[index].Pt()>6.)) break; | |
568 | } | |
f3050824 | 569 | } |
623ed0a0 | 570 | |
571 | //Set new Pt & Fill histogram accordingly | |
572 | maxPtJet1 = random * jetVectnew[index].Pt(); | |
573 | ||
574 | ||
575 | fhEleadingPt->Fill( maxPtJet1 ); | |
576 | ||
577 | if (fUseAliStack){//Try Stack Information to perform UE analysis | |
578 | ||
579 | AliStack* mcStack = mcEvent->Stack();//Load Stack | |
580 | Int_t nTracksMC = mcStack->GetNtrack(); | |
581 | for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) { | |
582 | //Cuts | |
583 | if(!(mcStack->IsPhysicalPrimary(iTracks))) continue; | |
584 | ||
585 | TParticle* mctrk = mcStack->Particle(iTracks); | |
586 | ||
587 | Double_t charge = mctrk->GetPDG()->Charge(); | |
588 | if (charge == 0) continue; | |
589 | ||
590 | if ( fUseSingleCharge ) { // Charge selection | |
591 | if ( fUsePositiveCharge && charge < 0.) continue; // keep Positives | |
592 | if ( !fUsePositiveCharge && charge > 0.) continue; // keep Negatives | |
593 | } | |
594 | ||
595 | //Kinematics cuts on particle | |
596 | if ((mctrk->Pt() < fTrackPtCut) || (TMath::Abs(mctrk->Eta()) > fTrackEtaCut )) continue; | |
597 | ||
598 | Bool_t isHadron = TMath::Abs(mctrk->GetPdgCode())==211 || | |
599 | TMath::Abs(mctrk->GetPdgCode())==2212 || | |
600 | TMath::Abs(mctrk->GetPdgCode())==321; | |
601 | ||
602 | if (!isHadron) continue; | |
603 | ||
604 | TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz()); | |
605 | ||
606 | Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad; | |
607 | if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi(); | |
608 | fhdNdEtaPhiDist->Fill( deltaPhi ); | |
609 | ||
610 | fhFullRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 ); | |
611 | ||
612 | //We are not interested on stack organization but don't loose track of info | |
613 | TVector3 tempVector = jetVectnew[0]; | |
614 | jetVectnew[0] = jetVectnew[index]; | |
615 | jetVectnew[index] = tempVector; | |
616 | ||
617 | Int_t region = IsTrackInsideRegion( jetVectnew, &partVect ); | |
618 | ||
619 | if (region > 0) { | |
620 | if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt(); | |
621 | sumPtRegionPosit += mctrk->Pt(); | |
622 | nTrackRegionPosit++; | |
623 | fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 ); | |
624 | } | |
625 | if (region < 0) { | |
626 | if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt(); | |
627 | sumPtRegionNegat += mctrk->Pt(); | |
628 | nTrackRegionNegat++; | |
629 | fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 ); | |
630 | } | |
631 | } // end loop stack Particles | |
632 | ||
633 | }else{//Try mc Particle | |
634 | ||
635 | TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles"); | |
636 | ||
637 | Int_t ntrks = farray->GetEntries(); | |
638 | if (fDebug>1) AliInfo(Form("In UE MC analysis tracks %d \n",ntrks)); | |
639 | for(Int_t i =0 ; i < ntrks; i++){ | |
640 | AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i); | |
641 | //Cuts | |
642 | if (!(mctrk->IsPhysicalPrimary())) continue; | |
643 | //if (!(mctrk->IsPrimary())) continue; | |
644 | ||
645 | if (mctrk->Charge() == 0 || mctrk->Charge()==-99) continue; | |
646 | ||
647 | if (mctrk->Pt() < fTrackPtCut ) continue; | |
648 | if( TMath::Abs(mctrk->Eta()) > fTrackEtaCut ) continue; | |
649 | ||
650 | Bool_t isHadron = TMath::Abs(mctrk->GetPdgCode())==211 || | |
651 | TMath::Abs(mctrk->GetPdgCode())==2212 || | |
652 | TMath::Abs(mctrk->GetPdgCode())==321; | |
653 | ||
654 | if (!isHadron) continue; | |
655 | ||
656 | TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz()); | |
657 | ||
658 | Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad; | |
659 | if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi(); | |
660 | fhdNdEtaPhiDist->Fill( deltaPhi ); | |
661 | ||
662 | fhFullRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 ); | |
663 | ||
664 | //We are not interested on stack organization but don't loose track of info | |
665 | TVector3 tempVector = jetVectnew[0]; | |
666 | jetVectnew[0] = jetVectnew[index]; | |
667 | jetVectnew[index] = tempVector; | |
668 | ||
669 | Int_t region = IsTrackInsideRegion( jetVectnew, &partVect ); | |
670 | ||
671 | if (region > 0) { | |
672 | if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt(); | |
673 | sumPtRegionPosit += mctrk->Pt(); | |
674 | nTrackRegionPosit++; | |
675 | fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 ); | |
676 | } | |
677 | if (region < 0) { | |
678 | if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt(); | |
679 | sumPtRegionNegat += mctrk->Pt(); | |
680 | nTrackRegionNegat++; | |
681 | fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 ); | |
682 | } | |
683 | ||
684 | }//end loop AliAODMCParticle tracks | |
f3050824 | 685 | } |
686 | } | |
f3050824 | 687 | //How quantities will be sorted before Fill Min and Max Histogram |
688 | // 1=Plots will be CDF-like | |
689 | // 2=Plots will be Marchesini-like | |
6ef5bfa4 | 690 | if( fOrdering == 1 ) { |
691 | if( sumPtRegionPosit > sumPtRegionNegat ) { | |
692 | FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg ); | |
f3050824 | 693 | } else { |
6ef5bfa4 | 694 | FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg ); |
f3050824 | 695 | } |
696 | if (nTrackRegionPosit > nTrackRegionNegat ) { | |
6ef5bfa4 | 697 | FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg ); |
f3050824 | 698 | } else { |
6ef5bfa4 | 699 | FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg ); |
f3050824 | 700 | } |
6ef5bfa4 | 701 | } else if( fOrdering == 2 ) { |
f3050824 | 702 | if (sumPtRegionPosit > sumPtRegionNegat) { |
6ef5bfa4 | 703 | FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg ); |
704 | FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg ); | |
f3050824 | 705 | } else { |
6ef5bfa4 | 706 | FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg ); |
707 | FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg ); | |
f3050824 | 708 | } |
709 | } | |
710 | ||
711 | Double_t avePosRegion = (nTrackRegionPosit) ? sumPtRegionPosit/nTrackRegionPosit : 0.; | |
712 | Double_t aveNegRegion = (nTrackRegionNegat) ? sumPtRegionNegat/nTrackRegionNegat : 0.; | |
6ef5bfa4 | 713 | if( avePosRegion > aveNegRegion ) { |
714 | FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg ); | |
f3050824 | 715 | } else { |
6ef5bfa4 | 716 | FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg ); |
f3050824 | 717 | } |
6ef5bfa4 | 718 | |
f3050824 | 719 | fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion ); |
6ef5bfa4 | 720 | |
721 | // Compute pedestal like magnitudes | |
722 | fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/(2.0*fAreaReg)); | |
723 | fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/(2.0*fAreaReg)); | |
724 | ||
f3050824 | 725 | } |
726 | ||
f3050824 | 727 | //____________________________________________________________________ |
728 | void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ) | |
729 | { | |
6ef5bfa4 | 730 | // Fill sumPt of control regions |
731 | ||
732 | // Max cone | |
733 | fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax ); | |
734 | // Min cone | |
735 | fhRegionSumPtMinVsEt->Fill( leadingE, ptMin ); | |
736 | // MAke distributions for UE comparison with MB data | |
737 | fhMinRegSumPt->Fill(ptMin); | |
f3050824 | 738 | } |
739 | ||
740 | //____________________________________________________________________ | |
741 | void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ) | |
742 | { | |
6ef5bfa4 | 743 | // Fill average particle Pt of control regions |
744 | ||
745 | // Max cone | |
746 | fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax ); | |
747 | // Min cone | |
748 | fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin ); | |
749 | // MAke distributions for UE comparison with MB data | |
750 | fhMinRegAvePt->Fill(ptMin); | |
f3050824 | 751 | } |
752 | ||
753 | //____________________________________________________________________ | |
754 | void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin ) | |
755 | { | |
6ef5bfa4 | 756 | // Fill Nch multiplicity of control regions |
757 | ||
758 | // Max cone | |
759 | fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax ); | |
760 | fhRegionMultMax->Fill( nTrackPtmax ); | |
761 | // Min cone | |
762 | fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin ); | |
763 | fhRegionMultMin->Fill( nTrackPtmin ); | |
764 | // MAke distributions for UE comparison with MB data | |
765 | fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin); | |
f3050824 | 766 | } |
767 | ||
768 | //____________________________________________________________________ | |
769 | Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect) | |
770 | { | |
771 | // return de region in delta phi | |
772 | // -1 negative delta phi | |
773 | // 1 positive delta phi | |
774 | // 0 outside region | |
775 | static const Double_t k60rad = 60.*TMath::Pi()/180.; | |
776 | static const Double_t k120rad = 120.*TMath::Pi()/180.; | |
777 | ||
778 | Int_t region = 0; | |
779 | if( fRegionType == 1 ) { | |
780 | if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0; | |
781 | // transverse regions | |
782 | if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; | |
783 | if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1; | |
6ef5bfa4 | 784 | |
f3050824 | 785 | } else if( fRegionType == 2 ) { |
786 | // Cone regions | |
787 | Double_t deltaR = 0.; | |
788 | ||
789 | TVector3 positVect,negatVect; | |
623ed0a0 | 790 | if (fConePosition==1){ |
791 | positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2()); | |
792 | negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2()); | |
793 | }else if (fConePosition==2){ | |
794 | if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config"); | |
795 | positVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()+TMath::PiOver2()); | |
796 | negatVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()-TMath::PiOver2()); | |
797 | }else if (fConePosition==3){ | |
798 | if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config"); | |
799 | Double_t weightEta = jetVect[0].Eta() * jetVect[0].Pt()/(jetVect[0].Pt() + jetVect[1].Pt()) + | |
800 | jetVect[1].Eta() * jetVect[1].Pt()/(jetVect[0].Pt() + jetVect[1].Pt()); | |
801 | //Double_t weightEta = jetVect[0].Eta() * jetVect[0].Mag()/(jetVect[0].Mag() + jetVect[1].Mag()) + | |
802 | // jetVect[1].Eta() * jetVect[1].Mag()/(jetVect[0].Mag() + jetVect[1].Mag()); | |
803 | positVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()+TMath::PiOver2()); | |
804 | negatVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()-TMath::PiOver2()); | |
805 | } | |
f3050824 | 806 | if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) { |
6ef5bfa4 | 807 | region = 1; |
808 | deltaR = positVect.DrEtaPhi(*partVect); | |
f3050824 | 809 | } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) { |
6ef5bfa4 | 810 | region = -1; |
811 | deltaR = negatVect.DrEtaPhi(*partVect); | |
f3050824 | 812 | } |
813 | ||
814 | if (deltaR > fConeRadius) region = 0; | |
6ef5bfa4 | 815 | |
f3050824 | 816 | } else |
817 | AliError("Unknow region type"); | |
6ef5bfa4 | 818 | |
f3050824 | 819 | // For debug (to be removed) |
820 | //if( region != 0 ) fhValidRegion->Fill( partVect->Eta()-jetVect[0].Eta(), jetVect[0].DeltaPhi(*partVect) ); | |
821 | ||
822 | return region; | |
823 | } | |
824 | ||
6ef5bfa4 | 825 | |
826 | //____________________________________________________________________ | |
827 | TObjArray* AliAnalysisTaskUE::FindChargedParticleJets() | |
828 | { | |
829 | // Return a TObjArray of "charged particle jets" | |
830 | // | |
831 | // Charged particle jet definition from reference: | |
832 | // | |
833 | // "Charged jet evolution and the underlying event | |
834 | // in proton-antiproton collisions at 1.8 TeV" | |
835 | // PHYSICAL REVIEW D 65 092002, CDF Collaboration | |
836 | // | |
837 | // We define "jets" as circular regions in eta-phi space with | |
838 | // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ). | |
839 | // Our jet algorithm is as follows: | |
840 | // 1- Order all charged particles according to their pT . | |
841 | // 2- Start with the highest pT particle and include in the jet all | |
842 | // particles within the radius R=0.7 considering each particle | |
843 | // in the order of decreasing pT and recalculating the centroid | |
844 | // of the jet after each new particle is added to the jet . | |
845 | // 3- Go to the next highest pT particle not already included in | |
846 | // a jet and add to the jet all particles not already included in | |
847 | // a jet within R=0.7. | |
848 | // 4- Continue until all particles are in a jet. | |
849 | // We define the transverse momentum of the jet to be | |
850 | // the scalar pT sum of all the particles within the jet, where pT | |
851 | // is measured with respect to the beam axis | |
852 | ||
853 | // 1 - Order all charged particles according to their pT . | |
854 | Int_t nTracks = fAOD->GetNTracks(); | |
855 | if( !nTracks ) return 0; | |
856 | TObjArray tracks(nTracks); | |
857 | ||
858 | for (Int_t ipart=0; ipart<nTracks; ++ipart) { | |
859 | AliAODTrack* part = fAOD->GetTrack( ipart ); | |
860 | if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection | |
861 | if( !part->Charge() ) continue; | |
862 | if( part->Pt() < fTrackPtCut ) continue; | |
863 | tracks.AddLast(part); | |
864 | } | |
865 | QSortTracks( tracks, 0, tracks.GetEntriesFast() ); | |
866 | ||
867 | nTracks = tracks.GetEntriesFast(); | |
868 | if( !nTracks ) return 0; | |
869 | TObjArray *jets = new TObjArray(nTracks); | |
870 | TIter itrack(&tracks); | |
871 | while( nTracks ) { | |
872 | // 2- Start with the highest pT particle ... | |
873 | Float_t px,py,pz,pt; | |
874 | AliAODTrack* track = (AliAODTrack*)itrack.Next(); | |
875 | if( !track ) continue; | |
876 | px = track->Px(); | |
877 | py = track->Py(); | |
878 | pz = track->Pz(); | |
879 | pt = track->Pt(); // Use the energy member to store Pt | |
880 | jets->AddLast( new TLorentzVector(px, py, pz, pt) ); | |
881 | tracks.Remove( track ); | |
882 | TLorentzVector* jet = (TLorentzVector*)jets->Last(); | |
883 | // 3- Go to the next highest pT particle not already included... | |
884 | AliAODTrack* track1; | |
885 | while ( (track1 = (AliAODTrack*)(itrack.Next())) ) { | |
886 | Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-track1->Phi()); | |
887 | Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) + | |
888 | dphi*dphi ); | |
889 | if( r < fConeRadius ) { | |
de6f8090 | 890 | Double_t fPt = jet->E()+track1->Pt(); // Scalar sum of Pt |
6ef5bfa4 | 891 | // recalculating the centroid |
623ed0a0 | 892 | Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt; |
893 | Double_t phi = jet->Phi()*jet->E()/fPt + track1->Phi()*track1->Pt()/fPt; | |
de6f8090 | 894 | jet->SetPtEtaPhiE( 1., eta, phi, fPt ); |
6ef5bfa4 | 895 | tracks.Remove( track1 ); |
896 | } | |
897 | } | |
898 | ||
899 | tracks.Compress(); | |
900 | nTracks = tracks.GetEntries(); | |
901 | // 4- Continue until all particles are in a jet. | |
902 | itrack.Reset(); | |
903 | } // end while nTracks | |
904 | ||
905 | // Convert to AODjets.... | |
906 | Int_t njets = jets->GetEntriesFast(); | |
907 | TObjArray* aodjets = new TObjArray(njets); | |
908 | aodjets->SetOwner(kTRUE); | |
909 | for(Int_t ijet=0; ijet<njets; ++ijet) { | |
910 | TLorentzVector* jet = (TLorentzVector*)jets->At(ijet); | |
623ed0a0 | 911 | if (jet->E() < fChJetPtMin) continue; |
6ef5bfa4 | 912 | Float_t px, py,pz,en; // convert to 4-vector |
913 | px = jet->E() * TMath::Cos(jet->Phi()); // Pt * cos(phi) | |
914 | py = jet->E() * TMath::Sin(jet->Phi()); // Pt * sin(phi) | |
915 | pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta()))); | |
916 | en = TMath::Sqrt(px * px + py * py + pz * pz); | |
917 | aodjets->AddLast( new AliAODJet(px, py, pz, en) ); | |
918 | } | |
919 | // Order jets according to their pT . | |
920 | QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() ); | |
921 | ||
922 | // debug | |
923 | if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets)); | |
924 | ||
925 | return aodjets; | |
926 | } | |
927 | ||
928 | //____________________________________________________________________ | |
929 | void AliAnalysisTaskUE::QSortTracks(TObjArray &a, Int_t first, Int_t last) | |
930 | { | |
931 | // Sort array of TObjArray of tracks by Pt using a quicksort algorithm. | |
932 | ||
933 | static TObject *tmp; | |
934 | static int i; // "static" to save stack space | |
935 | int j; | |
936 | ||
937 | while (last - first > 1) { | |
938 | i = first; | |
939 | j = last; | |
940 | for (;;) { | |
941 | while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() ) | |
942 | ; | |
943 | while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() ) | |
944 | ; | |
945 | if (i >= j) | |
946 | break; | |
947 | ||
948 | tmp = a[i]; | |
949 | a[i] = a[j]; | |
950 | a[j] = tmp; | |
951 | } | |
952 | if (j == first) { | |
953 | ++first; | |
954 | continue; | |
955 | } | |
956 | tmp = a[first]; | |
957 | a[first] = a[j]; | |
958 | a[j] = tmp; | |
959 | if (j - first < last - (j + 1)) { | |
960 | QSortTracks(a, first, j); | |
961 | first = j + 1; // QSortTracks(j + 1, last); | |
962 | } else { | |
963 | QSortTracks(a, j + 1, last); | |
964 | last = j; // QSortTracks(first, j); | |
965 | } | |
966 | } | |
967 | } | |
968 | ||
969 | //____________________________________________________________________ | |
970 | void AliAnalysisTaskUE::SetRegionArea(TVector3 *jetVect) | |
971 | { | |
972 | Double_t fAreaCorrFactor=0.; | |
973 | Double_t deltaEta = 0.; | |
974 | if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.; | |
975 | else if (fRegionType==2){ | |
976 | deltaEta = 0.9-TMath::Abs(jetVect[0].Eta()); | |
977 | if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius; | |
978 | else{ | |
979 | fAreaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) - | |
980 | (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius)); | |
981 | fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-fAreaCorrFactor; | |
982 | } | |
983 | }else AliWarning("Unknown Rgion Type"); | |
623ed0a0 | 984 | if (fDebug>10) AliInfo(Form("\n dEta=%5.3f Angle =%5.3f Region Area = %5.3f Corr Factor=%5.4f \n",deltaEta,TMath::ACos(deltaEta/fConeRadius),fAreaReg,fAreaCorrFactor)); |
6ef5bfa4 | 985 | } |
986 | ||
f3050824 | 987 | //____________________________________________________________________ |
988 | void AliAnalysisTaskUE::CreateHistos() | |
989 | { | |
990 | fListOfHistos = new TList(); | |
991 | ||
992 | ||
993 | fhNJets = new TH1F("fhNJets", "n Jet", 10, 0, 10); | |
994 | fhNJets->SetXTitle("# of jets"); | |
995 | fhNJets->Sumw2(); | |
f3050824 | 996 | fListOfHistos->Add( fhNJets ); // At(0) |
997 | ||
998 | fhEleadingPt = new TH1F("hEleadingPt", "Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); | |
999 | fhEleadingPt->SetXTitle("P_{T} (GeV/c)"); | |
1000 | fhEleadingPt->SetYTitle("dN/dP_{T}"); | |
1001 | fhEleadingPt->Sumw2(); | |
1002 | fListOfHistos->Add( fhEleadingPt ); // At(1) | |
1003 | ||
1004 | fhMinRegPtDist = new TH1F("hMinRegPtDist", "P_{T} distribution in Min zone", 50,0.,20.); | |
1005 | fhMinRegPtDist->SetXTitle("P_{T} (GeV/c)"); | |
1006 | fhMinRegPtDist->SetYTitle("dN/dP_{T}"); | |
1007 | fhMinRegPtDist->Sumw2(); | |
1008 | fListOfHistos->Add( fhMinRegPtDist ); // At(2) | |
1009 | ||
1010 | fhRegionMultMin = new TH1F("hRegionMultMin", "N_{ch}^{90, min}", 21, -0.5, 20.5); | |
1011 | fhRegionMultMin->SetXTitle("N_{ch tracks}"); | |
1012 | fhRegionMultMin->Sumw2(); | |
1013 | fListOfHistos->Add( fhRegionMultMin ); // At(3) | |
1014 | ||
1015 | fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT", 50, 0., 20.); | |
1016 | fhMinRegAvePt->SetXTitle("P_{T} (GeV/c)"); | |
1017 | fhMinRegAvePt->Sumw2(); | |
1018 | fListOfHistos->Add( fhMinRegAvePt ); // At(4) | |
1019 | ||
1020 | fhMinRegSumPt = new TH1F("hMinRegSumPt", "#Sigma p_{T} ", 50, 0., 20.); | |
1021 | fhMinRegSumPt->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})"); | |
1022 | fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)"); | |
1023 | fhMinRegSumPt->Sumw2(); | |
1024 | fListOfHistos->Add( fhMinRegSumPt ); // At(5) | |
6ef5bfa4 | 1025 | |
f3050824 | 1026 | fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ", 50, 0., 20.); |
1027 | fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})"); | |
1028 | fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)"); | |
1029 | fhMinRegMaxPtPart->Sumw2(); | |
1030 | fListOfHistos->Add( fhMinRegMaxPtPart ); // At(6) | |
1031 | ||
1032 | fhMinRegSumPtvsMult = new TH1F("hMinRegSumPtvsMult", "#Sigma p_{T} vs. Multiplicity ", 21, -0.5, 20.5); | |
1033 | fhMinRegSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)"); | |
1034 | fhMinRegSumPtvsMult->SetXTitle("N_{charge}"); | |
1035 | fhMinRegSumPtvsMult->Sumw2(); | |
1036 | fListOfHistos->Add( fhMinRegSumPtvsMult ); // At(7); | |
6ef5bfa4 | 1037 | |
de6f8090 | 1038 | fhdNdEtaPhiDist = new TH1F("hdNdEtaPhiDist", "Charge particle density |#eta|< 1 vs #Delta#phi", 120, 0., 2.*TMath::Pi()); |
1039 | fhdNdEtaPhiDist->SetXTitle("#Delta#phi"); | |
1040 | fhdNdEtaPhiDist->SetYTitle("dN_{ch}/d#etad#phi"); | |
1041 | fhdNdEtaPhiDist->Sumw2(); | |
1042 | fListOfHistos->Add( fhdNdEtaPhiDist ); // At(8) | |
6ef5bfa4 | 1043 | |
f3050824 | 1044 | // Can be use to get part pt distribution for differente Jet Pt bins |
1045 | fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", "dN/dP_{T} |#eta|<1 vs Leading Jet P_{T}", | |
6ef5bfa4 | 1046 | 50,0.,50., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); |
f3050824 | 1047 | fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}"); |
1048 | fhFullRegPartPtDistVsEt->SetXTitle("p_{T}"); | |
1049 | fhFullRegPartPtDistVsEt->Sumw2(); | |
1050 | fListOfHistos->Add( fhFullRegPartPtDistVsEt ); // At(9) | |
6ef5bfa4 | 1051 | |
1052 | // Can be use to get part pt distribution for differente Jet Pt bins | |
f3050824 | 1053 | fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", "dN/dP_{T} in tranvese regions |#eta|<1 vs Leading Jet P_{T}", |
6ef5bfa4 | 1054 | 50,0.,50., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); |
f3050824 | 1055 | fhTransRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}"); |
1056 | fhTransRegPartPtDistVsEt->SetXTitle("p_{T}"); | |
1057 | fhTransRegPartPtDistVsEt->Sumw2(); | |
1058 | fListOfHistos->Add( fhTransRegPartPtDistVsEt ); // At(10) | |
1059 | ||
1060 | ||
1061 | fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt", "P_{T}^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); | |
1062 | fhRegionSumPtMaxVsEt->SetXTitle("P_{T} (GeV/c)"); | |
1063 | fhRegionSumPtMaxVsEt->Sumw2(); | |
1064 | fListOfHistos->Add( fhRegionSumPtMaxVsEt ); // At(11) | |
6ef5bfa4 | 1065 | |
f3050824 | 1066 | fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt", "P_{T}^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); |
1067 | fhRegionSumPtMinVsEt->SetXTitle("P_{T} (GeV/c)"); | |
1068 | fhRegionSumPtMinVsEt->Sumw2(); | |
1069 | fListOfHistos->Add( fhRegionSumPtMinVsEt ); // At(12) | |
1070 | ||
1071 | fhRegionMultMax = new TH1I("hRegionMultMax", "N_{ch}^{90, max}", 21, -0.5, 20.5); | |
1072 | fhRegionMultMax->SetXTitle("N_{ch tracks}"); | |
1073 | fhRegionMultMax->Sumw2(); | |
1074 | fListOfHistos->Add( fhRegionMultMax ); // At(13) | |
6ef5bfa4 | 1075 | |
f3050824 | 1076 | fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt", "N_{ch}^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); |
1077 | fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)"); | |
1078 | fhRegionMultMaxVsEt->Sumw2(); | |
1079 | fListOfHistos->Add( fhRegionMultMaxVsEt ); // At(14) | |
6ef5bfa4 | 1080 | |
f3050824 | 1081 | fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt", "N_{ch}^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); |
1082 | fhRegionMultMinVsEt->SetXTitle("E (GeV/c)"); | |
1083 | fhRegionMultMinVsEt->Sumw2(); | |
1084 | fListOfHistos->Add( fhRegionMultMinVsEt ); // At(15) | |
6ef5bfa4 | 1085 | |
f3050824 | 1086 | fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); |
1087 | fhRegionAveSumPtVsEt->SetXTitle("P_{T} (GeV/c)"); | |
1088 | fhRegionAveSumPtVsEt->Sumw2(); | |
1089 | fListOfHistos->Add( fhRegionAveSumPtVsEt ); // At(16) | |
6ef5bfa4 | 1090 | |
f3050824 | 1091 | fhRegionDiffSumPtVsEt= new TH1F("hRegionPtDiffVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); |
1092 | fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)"); | |
1093 | fhRegionDiffSumPtVsEt->Sumw2(); | |
1094 | fListOfHistos->Add( fhRegionDiffSumPtVsEt ); // At(17) | |
1095 | ||
1096 | fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); | |
1097 | fhRegionAvePartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)"); | |
1098 | fhRegionAvePartPtMaxVsEt->Sumw2(); | |
1099 | fListOfHistos->Add( fhRegionAvePartPtMaxVsEt ); // At(18) | |
6ef5bfa4 | 1100 | |
f3050824 | 1101 | fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); |
1102 | fhRegionAvePartPtMinVsEt->SetXTitle("P_{T} (GeV/c)"); | |
1103 | fhRegionAvePartPtMinVsEt->Sumw2(); | |
1104 | fListOfHistos->Add( fhRegionAvePartPtMinVsEt ); // At(19) | |
6ef5bfa4 | 1105 | |
f3050824 | 1106 | fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); |
1107 | fhRegionMaxPartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)"); | |
1108 | fhRegionMaxPartPtMaxVsEt->Sumw2(); | |
1109 | fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt ); // At(20) | |
6ef5bfa4 | 1110 | |
623ed0a0 | 1111 | fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1); |
1112 | fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>"); | |
1113 | fh1Xsec->Sumw2(); | |
1114 | fListOfHistos->Add( fh1Xsec ); //At(21) | |
1115 | ||
1116 | fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1); | |
1117 | fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}"); | |
1118 | fh1Trials->Sumw2(); | |
1119 | fListOfHistos->Add( fh1Trials ); //At(22) | |
1e996f55 | 1120 | |
6ef5bfa4 | 1121 | fSettingsTree = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation"); |
1e996f55 | 1122 | fSettingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D"); |
1123 | fSettingsTree->Branch("fJet1EtaCut", &fJet1EtaCut, "LeadJetEtaCut/D"); | |
1124 | fSettingsTree->Branch("fJet2DeltaPhiCut", &fJet2DeltaPhiCut, "DeltaPhi/D"); | |
1125 | fSettingsTree->Branch("fJet2RatioPtCut", &fJet2RatioPtCut, "Jet2Ratio/D"); | |
1126 | fSettingsTree->Branch("fJet3PtCut", &fJet3PtCut, "Jet3PtCut/D"); | |
1127 | fSettingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D"); | |
1128 | fSettingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D"); | |
1129 | fSettingsTree->Branch("fAnaType", &fAnaType, "Ana/I"); | |
1130 | fSettingsTree->Branch("fRegionType", &fRegionType,"Reg/I"); | |
1131 | fSettingsTree->Branch("fOrdering", &fOrdering,"OrderMeth/I"); | |
1132 | fSettingsTree->Branch("fUseChPartJet", &fUseChPartJet,"UseChPart/O"); | |
1133 | fSettingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O"); | |
1134 | fSettingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O"); | |
1135 | fSettingsTree->Fill(); | |
1136 | ||
1137 | ||
623ed0a0 | 1138 | fListOfHistos->Add( fSettingsTree ); // At(23) |
6ef5bfa4 | 1139 | |
1140 | /* | |
1141 | // For debug region selection | |
1142 | fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi", | |
1143 | 20, -2.,2., 62, -TMath::Pi(), TMath::Pi()); | |
1144 | fhValidRegion->SetYTitle("#Delta#phi"); | |
1145 | fhValidRegion->Sumw2(); | |
1146 | fListOfHistos->Add( fhValidRegion ); // At(15) | |
1147 | */ | |
f3050824 | 1148 | } |
f3050824 | 1149 | |
6ef5bfa4 | 1150 | //____________________________________________________________________ |
1151 | void AliAnalysisTaskUE::Terminate(Option_t */*option*/) | |
1152 | { | |
1153 | // Terminate analysis | |
1154 | // | |
1155 | ||
1156 | //Save Analysis Settings | |
1e996f55 | 1157 | // WriteSettings(); |
6ef5bfa4 | 1158 | |
1159 | // Normalize histos to region area TODO: | |
1160 | // Normalization done at Analysis, taking into account | |
1161 | // area variations on per-event basis (cone case) | |
1162 | ||
1163 | // Draw some Test plot to the screen | |
1164 | if (!gROOT->IsBatch()) { | |
1165 | ||
1166 | // Update pointers reading them from the output slot | |
1167 | fListOfHistos = dynamic_cast<TList*> (GetOutputData(0)); | |
1168 | if( !fListOfHistos ) { | |
1169 | AliError("Histogram List is not available"); | |
1170 | return; | |
1171 | } | |
1172 | fhEleadingPt = (TH1F*)fListOfHistos->At(1); | |
de6f8090 | 1173 | fhdNdEtaPhiDist = (TH1F*)fListOfHistos->At(8); |
6ef5bfa4 | 1174 | fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11); |
1175 | fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12); | |
1176 | fhRegionMultMaxVsEt = (TH1F*)fListOfHistos->At(14); | |
1177 | fhRegionMultMinVsEt = (TH1F*)fListOfHistos->At(15); | |
1178 | fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16); | |
1179 | ||
1180 | fhNJets = (TH1F*)fListOfHistos->At(0); | |
1181 | ||
623ed0a0 | 1182 | //fhValidRegion = (TH2F*)fListOfHistos->At(21); |
1183 | ||
6ef5bfa4 | 1184 | // Canvas |
623ed0a0 | 1185 | TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1100,700); |
6ef5bfa4 | 1186 | c1->Divide(2,2); |
1187 | c1->cd(1); | |
1188 | TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); | |
1189 | h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1); | |
1190 | //h1r->Scale( areafactor ); | |
1191 | h1r->SetMarkerStyle(20); | |
1192 | h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); | |
1193 | h1r->SetYTitle("P_{T}^{90, max}"); | |
1194 | h1r->DrawCopy("p"); | |
1195 | ||
1196 | c1->cd(2); | |
1197 | ||
1198 | TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); | |
1199 | h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1); | |
1200 | //h2r->Scale( areafactor ); | |
1201 | h2r->SetMarkerStyle(20); | |
1202 | h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); | |
1203 | h2r->SetYTitle("P_{T}^{90, min}"); | |
1204 | h2r->DrawCopy("p"); | |
1205 | ||
1206 | c1->cd(3); | |
1207 | TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); | |
1208 | h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1); | |
1209 | h4r->Scale(2.); // make average | |
1210 | //h4r->Scale( areafactor ); | |
1211 | h4r->SetYTitle("#DeltaP_{T}^{90}"); | |
1212 | h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); | |
1213 | h4r->SetMarkerStyle(20); | |
1214 | h4r->DrawCopy("p"); | |
1215 | ||
1216 | c1->cd(4); | |
1217 | TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); | |
1218 | TH1F *h6r = new TH1F("hRegionMultMinVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); | |
1219 | ||
1220 | h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1); | |
1221 | h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1); | |
1222 | //h5r->Scale( areafactor ); | |
1223 | //h6r->Scale( areafactor ); | |
1224 | h5r->SetYTitle("N_{Tracks}^{90}"); | |
1225 | h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); | |
1226 | h5r->SetMarkerStyle(20); | |
1227 | h5r->DrawCopy("p"); | |
1228 | h6r->SetMarkerStyle(21); | |
1229 | h6r->SetMarkerColor(2); | |
1230 | h6r->SetYTitle("N_{Tracks}^{90}"); | |
1231 | h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); | |
1232 | h6r->DrawCopy("p same"); | |
1233 | c1->Update(); | |
1234 | ||
623ed0a0 | 1235 | //Get Normalization |
1236 | fh1Xsec = (TProfile*)fListOfHistos->At(21); | |
1237 | fh1Trials = (TH1F*)fListOfHistos->At(22); | |
1238 | ||
1239 | Double_t xsec = fh1Xsec->GetBinContent(1); | |
1240 | Double_t ntrials = fh1Trials->GetBinContent(1); | |
1241 | Double_t normFactor = xsec/ntrials; | |
1242 | Printf("xSec %f nTrials %f Norm %f \n",xsec,ntrials,normFactor); | |
1243 | ||
1244 | ||
6ef5bfa4 | 1245 | TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800); |
1246 | c2->Divide(2,2); | |
1247 | c2->cd(1); | |
1248 | fhEleadingPt->SetMarkerStyle(20); | |
1249 | fhEleadingPt->SetMarkerColor(2); | |
623ed0a0 | 1250 | fhEleadingPt->Scale(normFactor); |
1251 | //fhEleadingPt->Draw("p"); | |
6ef5bfa4 | 1252 | fhEleadingPt->DrawCopy("p"); |
1253 | gPad->SetLogy(); | |
1254 | ||
1255 | c2->cd(2); | |
de6f8090 | 1256 | fhdNdEtaPhiDist->SetMarkerStyle(20); |
1257 | fhdNdEtaPhiDist->SetMarkerColor(2); | |
1258 | fhdNdEtaPhiDist->DrawCopy("p"); | |
6ef5bfa4 | 1259 | gPad->SetLogy(); |
1260 | ||
1261 | c2->cd(3); | |
1262 | fhNJets->DrawCopy(); | |
1263 | ||
623ed0a0 | 1264 | //c2->cd(4); |
1265 | //fhValidRegion->DrawCopy("p"); | |
1266 | ||
6ef5bfa4 | 1267 | //fhTransRegPartPtDist = (TH1F*)fListOfHistos->At(2); |
1268 | fhRegionMultMin = (TH1F*)fListOfHistos->At(3); | |
1269 | fhMinRegAvePt = (TH1F*)fListOfHistos->At(4); | |
1270 | fhMinRegSumPt = (TH1F*)fListOfHistos->At(5); | |
1271 | //fhMinRegMaxPtPart = (TH1F*)fListOfHistos->At(6); | |
1272 | fhMinRegSumPtvsMult = (TH1F*)fListOfHistos->At(7); | |
1273 | ||
1274 | // Canvas | |
623ed0a0 | 1275 | TCanvas* c3 = new TCanvas("c3"," pT dist",160,160,1200,800); |
6ef5bfa4 | 1276 | c3->Divide(2,2); |
1277 | c3->cd(1); | |
1278 | /*fhTransRegPartPtDist->SetMarkerStyle(20); | |
1279 | fhTransRegPartPtDist->SetMarkerColor(2); | |
1280 | fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries()); | |
1281 | fhTransRegPartPtDist->DrawCopy("p"); | |
1282 | gPad->SetLogy(); | |
1283 | */ | |
1284 | c3->cd(2); | |
1285 | fhMinRegSumPt->SetMarkerStyle(20); | |
1286 | fhMinRegSumPt->SetMarkerColor(2); | |
1287 | //fhMinRegSumPt->Scale(areafactor); | |
1288 | fhMinRegSumPt->DrawCopy("p"); | |
1289 | gPad->SetLogy(); | |
1290 | ||
1291 | c3->cd(3); | |
1292 | fhMinRegAvePt->SetMarkerStyle(20); | |
1293 | fhMinRegAvePt->SetMarkerColor(2); | |
1294 | //fhMinRegAvePt->Scale(areafactor); | |
1295 | fhMinRegAvePt->DrawCopy("p"); | |
1296 | gPad->SetLogy(); | |
1297 | ||
1298 | c3->cd(4); | |
1299 | ||
1300 | TH1F *h7r = new TH1F("hRegionMultMinVsMult", "", 21, -0.5, 20.5); | |
1301 | ||
1302 | h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1); | |
1303 | ||
1304 | h7r->SetMarkerStyle(20); | |
1305 | h7r->SetMarkerColor(2); | |
1306 | h7r->DrawCopy("p"); | |
1307 | ||
1308 | c3->Update(); | |
1309 | ||
1310 | ||
1311 | /* c2->cd(3); | |
1312 | fhValidRegion->SetMarkerStyle(20); | |
1313 | fhValidRegion->SetMarkerColor(2); | |
1314 | fhValidRegion->DrawCopy("p"); | |
1315 | */ | |
1316 | c2->Update(); | |
1317 | } else { | |
1318 | AliInfo(" Batch mode, not histograms will be shown..."); | |
1319 | } | |
1320 | ||
1321 | if( fDebug > 1 ) | |
1322 | AliInfo("End analysis"); | |
1323 | ||
1324 | } | |
f3050824 | 1325 | |
6ef5bfa4 | 1326 | void AliAnalysisTaskUE::WriteSettings() |
1e996f55 | 1327 | { |
6ef5bfa4 | 1328 | if (fDebug>5){ |
1329 | AliInfo(" All Analysis Settings in Saved Tree"); | |
1330 | fSettingsTree->Scan(); | |
1331 | } | |
1332 | } |