1-New correlation analysis for particle and jet-JETAN correlation; 2-Change GetEntrie...
[u/mrichter/AliRoot.git] / PWG4 / AliAnalysisTaskJetSpectrum.cxx
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
52 ClassImp(AliAnalysisTaskJetSpectrum)
53
54 AliAnalysisTaskJetSpectrum::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
83 AliAnalysisTaskJetSpectrum::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
118 void 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
291 void 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
302 void 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
487 void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
488 {
489 // Terminate analysis
490 //
491     if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
492 }
493
494
495 void 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