]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/AliAnalysisTaskJetSpectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add Spectrum Task and helper class, some minor mods
[u/mrichter/AliRoot.git] / PWG4 / AliAnalysisTaskJetSpectrum.cxx
CommitLineData
cfff6259 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
17
18#include <TROOT.h>
19#include <TSystem.h>
20#include <TInterpreter.h>
21#include <TChain.h>
22#include <TFile.h>
23#include <TH1F.h>
24#include <TH2F.h>
25#include <TH3F.h>
26#include <TList.h>
27#include <TLorentzVector.h>
28#include <TClonesArray.h>
29#include "TDatabasePDG.h"
30
31#include "AliAnalysisTaskJetSpectrum.h"
32#include "AliAnalysisManager.h"
33#include "AliJetFinder.h"
34#include "AliJetReader.h"
35#include "AliJetReaderHeader.h"
36#include "AliUA1JetHeaderV1.h"
37#include "AliJet.h"
38#include "AliESDEvent.h"
39#include "AliAODEvent.h"
40#include "AliAODHandler.h"
41#include "AliAODTrack.h"
42#include "AliAODJet.h"
43#include "AliMCEventHandler.h"
44#include "AliMCEvent.h"
45#include "AliStack.h"
46#include "AliGenPythiaEventHeader.h"
47#include "AliJetKineReaderHeader.h"
48#include "AliGenCocktailEventHeader.h"
49
50#include "AliAnalysisHelperJetTasks.h"
51
52ClassImp(AliAnalysisTaskJetSpectrum)
53
54AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
55 fJetFinderRec(0x0),
56 fJetFinderGen(0x0),
57 fAOD(0x0),
58 fBranchRec("jets"),
59 fConfigRec("ConfigJets.C"),
60 fBranchGen(""),
61 fConfigGen(""),
62 fUseAODInput(kFALSE),
63 fUseExternalWeightOnly(kFALSE),
64 fAnalysisType(0),
65 fExternalWeight(1),
66 fh1PtHard(0x0),
67 fh1PtHard_NoW(0x0),
68 fh1PtHard_Trials(0x0),
69 fh1PtHard_Trials_NoW(0x0),
70 fh1NGenJets(0x0),
71 fh1NRecJets(0x0),
72 fHistList(0x0)
73{
74 // Default constructor
75 for(int ij = 0;ij<kMaxJets;++ij){
76 fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
77 fh2PtFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = 0;
78 fh3PtRecGenHard[ij] = fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3RecEtaPhiPt_NoFound[ij] = fh3MCEtaPhiPt[ij] = 0;
79 }
80
81}
82
83AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
84 AliAnalysisTaskSE(name),
85 fJetFinderRec(0x0),
86 fJetFinderGen(0x0),
87 fAOD(0x0),
88 fBranchRec("jets"),
89 fConfigRec("ConfigJets.C"),
90 fBranchGen(""),
91 fConfigGen(""),
92 fUseAODInput(kFALSE),
93 fUseExternalWeightOnly(kFALSE),
94 fAnalysisType(0),
95 fExternalWeight(1),
96 fh1PtHard(0x0),
97 fh1PtHard_NoW(0x0),
98 fh1PtHard_Trials(0x0),
99 fh1PtHard_Trials_NoW(0x0),
100 fh1NGenJets(0x0),
101 fh1NRecJets(0x0),
102 fHistList(0x0)
103{
104 // Default constructor
105 for(int ij = 0;ij<kMaxJets;++ij){
106 fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
107 fh2PtFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = 0;
108
109 fh3PtRecGenHard[ij] = fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3RecEtaPhiPt_NoFound[ij] = fh3MCEtaPhiPt[ij] = 0;
110 }
111
112 DefineOutput(1, TList::Class());
113}
114
115
116
117
118void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
119{
120
121 //
122 // Create the output container
123 //
124
125
126 // Connect the AOD
127
128 if(fUseAODInput){
129 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
130 if(!fAOD){
131 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
132 return;
133 }
134 }
135 else{
136 // assume that the AOD is in the general output...
137 fAOD = AODEvent();
138 if(!fAOD){
139 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
140 return;
141 }
142 }
143
144
145
146 if (fDebug > 1) printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
147
148 OpenFile(1);
149 if(!fHistList)fHistList = new TList();
150
151 Bool_t oldStatus = TH1::AddDirectoryStatus();
152 TH1::AddDirectory(kFALSE);
153
154 //
155 // Histogram
156
157 const Int_t nBinPt = 100;
158 Double_t binLimitsPt[nBinPt+1];
159 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
160 if(iPt == 0){
161 binLimitsPt[iPt] = 0.0;
162 }
163 else {// 1.0
164 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2;
165 }
166 }
167 const Int_t nBinEta = 22;
168 Double_t binLimitsEta[nBinEta+1];
169 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
170 if(iEta==0){
171 binLimitsEta[iEta] = -1.1;
172 }
173 else{
174 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
175 }
176 }
177
178 const Int_t nBinPhi = 360;
179 Double_t binLimitsPhi[nBinPhi+1];
180 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
181 if(iPhi==0){
182 binLimitsPhi[iPhi] = 0;
183 }
184 else{
185 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
186 }
187 }
188
189 const Int_t nBinFrag = 25;
190
191
192 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
193
194 fh1PtHard_NoW = new TH1F("fh1PtHard_NoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
195
196 fh1PtHard_Trials = new TH1F("fh1PtHard_Trials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
197
198 fh1PtHard_Trials_NoW = new TH1F("fh1PtHard_Trials_NoW","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
199
200 fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
201
202 fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
203
204
205 for(int ij = 0;ij<kMaxJets;++ij){
206 fh1E[ij] = new TH1F(Form("fh1E_j%d",ij),"Jet Energy;E_{jet} (GeV);N",nBinPt,binLimitsPt);
207 fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
208 fh1PtRecOut[ij] = new TH1F(Form("fh1PtRecOut_j%d",ij),"rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
209 fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
210 fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
211
212
213 fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
214 nBinPt,binLimitsPt,nBinPt,binLimitsPt);
215
216
217
218
219
220
221 fh3PtRecGenHard[ij] = new TH3F(Form("fh3PtRecGenHard_j%d",ij), "Pt hard vs. pt gen vs. pt rec;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
222
223
224
225 fh3PtRecGenHard_NoW[ij] = new TH3F(Form("fh3PtRecGenHard_NoW_j%d",ij), "Pt hard vs. pt gen vs. pt rec no weight;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
226
227
228 fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
229 nBinFrag,0.,1.,nBinPt,binLimitsPt);
230
231 fh2FragLn[ij] = new TH2F(Form("fh2FragLn_j%d",ij),"Jet Fragmentation Ln;#xi=ln(E_{jet}/E_{i});E_{jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
232 nBinFrag,0.,10.,nBinPt,binLimitsPt);
233
234 fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
235 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
236
237
238
239 fh3RecEtaPhiPt_NoGen[ij] = new TH3F(Form("fh3RecEtaPhiPt_NoGen_j%d",ij),"No generated for found jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
240 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
241
242
243 fh3RecEtaPhiPt_NoFound[ij] = new TH3F(Form("fh3RecEtaPhiPt_NoFound_g%d",ij),"No found for generated jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
244 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
245
246
247
248 fh3MCEtaPhiPt[ij] = new TH3F(Form("fh3MCEtaPhiPt_j%d",ij),"MC eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
249 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
250
251 }
252
253 const Int_t saveLevel = 1; // large save level more histos
254
255 if(saveLevel>0){
256 fHistList->Add(fh1PtHard);
257 fHistList->Add(fh1PtHard_NoW);
258 fHistList->Add(fh1PtHard_Trials);
259 fHistList->Add(fh1PtHard_Trials_NoW);
260 fHistList->Add(fh1NGenJets);
261 fHistList->Add(fh1NRecJets);
262 for(int ij = 0;ij<kMaxJets;++ij){
263 fHistList->Add(fh1E[ij]);
264 fHistList->Add(fh1PtRecIn[ij]);
265 fHistList->Add(fh1PtRecOut[ij]);
266 fHistList->Add(fh1PtGenIn[ij]);
267 fHistList->Add(fh1PtGenOut[ij]);
268 fHistList->Add(fh2PtFGen[ij]);
269 if(saveLevel>2){
270 fHistList->Add(fh3RecEtaPhiPt[ij]);
271 fHistList->Add(fh3RecEtaPhiPt_NoGen[ij]);
272 fHistList->Add(fh3RecEtaPhiPt_NoFound[ij]);
273 fHistList->Add(fh3MCEtaPhiPt[ij]);
274 }
275 }
276 }
277
278 // =========== Switch on Sumw2 for all histos ===========
279 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
280 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
281 if (h1){
282 // Printf("%s ",h1->GetName());
283 h1->Sumw2();
284 }
285 }
286
287 TH1::AddDirectory(oldStatus);
288
289}
290
291void AliAnalysisTaskJetSpectrum::Init()
292{
293 //
294 // Initialization
295 //
296
297 Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
298 if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
299
300}
301
302void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
303{
304 //
305 // Execute analysis for current event
306 //
307
308 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
309
310
311 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
312
313 if(!aodH){
314 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
315 return;
316 }
317
318
319 // aodH->SelectEvent(kTRUE);
320
321 // ========= These pointers need to be valid in any case =======
322
323
324 /*
325 AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
326 if(!jhRec){
327 Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
328 return;
329 }
330 */
331 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
332 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
333 if(!aodRecJets){
334 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
335 return;
336 }
337
338 // ==== General variables needed
339
340
341 // We use statice array, not to fragment the memory
342 AliAODJet genJets[kMaxJets];
343 Int_t nGenJets = 0;
344 AliAODJet recJets[kMaxJets];
345 Int_t nRecJets = 0;
346
347 Double_t eventW = 1;
348 Double_t ptHard = 0;
349 Double_t nTrials = 1; // Trials for MC trigger weigth for real data
350
351 if(fUseExternalWeightOnly){
352 eventW = fExternalWeight;
353 }
354
355
356 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
357 if((fAnalysisType&kAnaMC)==kAnaMC){
358 // this is the part we only use when we have MC information
359 AliMCEvent* mcEvent = MCEvent();
360 // AliStack *pStack = 0;
361 if(!mcEvent){
362 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
363 return;
364 }
365 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
366 if(!pythiaGenHeader){
367 return;
368 }
369
370 nTrials = pythiaGenHeader->Trials();
371 ptHard = pythiaGenHeader->GetPtHard();
372
373 if(!fUseExternalWeightOnly){
374 // case were we combine more than one p_T hard bin...
375 // eventW = AnalysisHelperCKB::GetXSectionWeight(pythiaGenHeader->GetPtHard(),fEnergy)*fExternalWeight;
376 }
377
378 // fetch the pythia generated jets only to be used here
379 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
380 AliAODJet pythiaGenJets[kMaxJets];
381
382 if(fBranchGen.Length()==0)nGenJets = nPythiaGenJets;
383 for(int ip = 0;ip < nPythiaGenJets;++ip){
384 if(ip>=kMaxJets)continue;
385 Float_t p[4];
386 pythiaGenHeader->TriggerJet(ip,p);
387 pythiaGenJets[ip].SetPxPyPzE(p[0],p[1],p[2],p[3]);
388 if(fBranchGen.Length()==0){
389 // if we have MC particles and we do not read from the aod branch
390 // use the pythia jets
391 genJets[ip].SetPxPyPzE(p[0],p[1],p[2],p[3]);
392 }
393 }
394
395
396 }// (fAnalysisType&kMC)==kMC)
397
398 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
399 fh1PtHard->Fill(ptHard,eventW);
400 fh1PtHard_NoW->Fill(ptHard,1);
401 fh1PtHard_Trials->Fill(ptHard,nTrials);
402
403
404 // If we set a second branch for the input jets fetch this
405 if(fBranchGen.Length()>0){
406 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
407 if(aodGenJets){
408 nGenJets = aodGenJets->GetEntries();
409 for(int ig = 0;ig < nGenJets;++ig){
410 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
411 if(!tmp)continue;
412 genJets[ig] = *tmp;
413 }
414 }
415 else{
416 Printf("%s:%d Generated jet branch %s not found",fBranchGen.Data());
417 }
418 }
419
420 fh1NGenJets->Fill(nGenJets);
421 // We do not want to exceed the maximum jet number
422 nGenJets = TMath::Min(nGenJets,kMaxJets);
423
424 // Fetch the reconstructed jets...
425
426
427 nRecJets = aodRecJets->GetEntries();
428 for(int ir = 0;ir < nRecJets;++ir){
429 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
430 if(!tmp)continue;
431 recJets[ir] = *tmp;
432 }
433
434 fh1NRecJets->Fill(nRecJets);
435 nRecJets = TMath::Min(nRecJets,kMaxJets);
436 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
437 // Relate the jets
438 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
439 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
440
441 for(int i = 0;i<kMaxJets;++i){
442 iGenIndex[i] = iRecIndex[i] = -1;
443 }
444
445
446 GetClosestJets(genJets,nGenJets,recJets,nRecJets,
447 iGenIndex,iRecIndex,fDebug);
448 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
449
450 // loop over reconstructed jets
451 for(int ir = 0;ir < nRecJets;++ir){
452 Double_t ptRec = recJets[ir].Pt();
453 Double_t phiRec = recJets[ir].Phi();
454 if(phiRec<0)phiRec+=TMath::Pi()*2.;
455 Double_t etaRec = recJets[ir].Eta();
456
457 fh1E[ir]->Fill(recJets[ir].E(),eventW);
458 fh1PtRecIn[ir]->Fill(ptRec,eventW);
459 fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
460 // Fill Correlation
461 Int_t ig = iGenIndex[ir];
462 if(ig>=0&&ig<nGenJets){
463 fh1PtRecOut[ir]->Fill(ptRec,eventW);
464 Double_t ptGen = genJets[ig].Pt();
465 fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
466 fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
467 fh3PtRecGenHard_NoW[ir]->Fill(ptRec,ptGen,ptHard,1);
468 }
469 }// loop over reconstructed jets
470
471
472 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
473 for(int ig = 0;ig < nGenJets;++ig){
474 Double_t ptGen = genJets[ig].Pt();
475 // Fill Correlation
476 fh1PtGenIn[ig]->Fill(ptGen,eventW);
477 Int_t ir = iRecIndex[ig];
478 if(ir>=0){
479 fh1PtGenOut[ig]->Fill(ptGen,eventW);
480 }
481 }// loop over reconstructed jets
482
483 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
484 PostData(1, fHistList);
485}
486
487void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
488{
489// Terminate analysis
490//
491 if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
492}
493
494
495void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,Int_t &nGenJets,
496 AliAODJet *recJets,Int_t &nRecJets,
497 Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
498
499 //
500 // Relate the two input jet Arrays
501 //
502
503 //
504 // The association has to be unique
505 // So check in two directions
506 // find the closest rec to a gen
507 // and check if there is no other rec which is closer
508 // Caveat: Close low energy/split jets may disturb this correlation
509
510 // Idea: search in two directions generated e.g (a--e) and rec (1--3)
511 // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
512 // in the end we have something like this
513 // 1 2 3
514 // ------------
515 // a| 3 2 0
516 // b| 0 1 0
517 // c| 0 0 3
518 // d| 0 0 1
519 // e| 0 0 1
520 // Topology
521 // 1 2
522 // a b
523 //
524 // d c
525 // 3 e
526 // Only entries with "3" match from both sides
527
528 Int_t iFlag[kMaxJets][kMaxJets];
529
530 for(int i = 0;i < kMaxJets;++i){
531 iRecIndex[i] = -1;
532 iGenIndex[i] = -1;
533 for(int j = 0;j < kMaxJets;++j)iFlag[i][j] = 0;
534 }
535
536 if(nRecJets==0)return;
537 if(nGenJets==0)return;
538
539 const Float_t maxDist = 1.4;
540 // find the closest distance
541 for(int ig = 0;ig<nGenJets;++ig){
542 Float_t dist = maxDist;
543 if(iDebug>1)Printf("Gen (%d) p_T %3.3f eta %3.3f ph %3.3f ",ig,genJets[ig].Pt(),genJets[ig].Eta(),genJets[ig].Phi());
544 for(int ir = 0;ir<nRecJets;++ir){
545 Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
546 if(iDebug>1)Printf("Rec (%d) p_T %3.3f eta %3.3f ph %3.3f ",ir,recJets[ir].Pt(),recJets[ir].Eta(),recJets[ir].Phi());
547 if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
548 if(dR<dist){
549 iRecIndex[ig] = ir;
550 dist = dR;
551 }
552 }
553 if(iRecIndex[ig]>=0)iFlag[ig][iRecIndex[ig]]+=1;
554 }
555 // other way around
556 for(int ir = 0;ir<nRecJets;++ir){
557 Float_t dist = maxDist;
558 for(int ig = 0;ig<nGenJets;++ig){
559 Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
560 if(dR<dist){
561 iGenIndex[ir] = ig;
562 dist = dR;
563 }
564 }
565 if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]][ir]+=2;
566 }
567
568 // check for "true" correlations
569
570 if(iDebug>1)Printf(">>>>>> Matrix");
571
572 for(int ig = 0;ig<nGenJets;++ig){
573 for(int ir = 0;ir<nRecJets;++ir){
574 // Print
575 if(iDebug>1)printf("%d ",iFlag[ig][ir]);
576 if(iFlag[ig][ir]==3){
577 iGenIndex[ir] = ig;
578 iRecIndex[ig] = ir;
579 }
580 else{
581 iGenIndex[ir] = iRecIndex[ig] = -1;
582 }
583 }
584 if(iDebug>1)printf("\n");
585 }
586}
587
588