Moved GetClosestJets() method to helper class
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisHelperJetTasks.cxx
CommitLineData
f3050824 1
2#include "TROOT.h"
3#include "TList.h"
4#include "AliMCEvent.h"
db6bcb0e 5#include "AliAODJet.h"
f3050824 6#include "AliStack.h"
7#include "AliGenEventHeader.h"
8#include "AliGenCocktailEventHeader.h"
9#include "AliGenPythiaEventHeader.h"
10#include <fstream>
11#include <iostream>
12#include "AliAnalysisHelperJetTasks.h"
13
14
15ClassImp(AliAnalysisHelperJetTasks)
16
17
18
19
20AliGenPythiaEventHeader* AliAnalysisHelperJetTasks::GetPythiaEventHeader(AliMCEvent *mcEvent){
21
22 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
23 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
24 if(!pythiaGenHeader){
25 // cocktail ??
26 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
27
28 if (!genCocktailHeader) {
29 Printf("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__);
30 return 0;
31 }
32 TList* headerList = genCocktailHeader->GetHeaders();
33 for (Int_t i=0; i<headerList->GetEntries(); i++) {
34 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
35 if (pythiaGenHeader)
36 break;
37 }
38 if(!pythiaGenHeader){
39 Printf("%s %d: PythiaHeader not found!",(char*)__FILE__,__LINE__);
40 return 0;
41 }
42 }
43 return pythiaGenHeader;
44
45}
46
47
48void AliAnalysisHelperJetTasks::PrintStack(AliMCEvent *mcEvent,Int_t iFirst,Int_t iLast,Int_t iMaxPrint){
49
50 AliStack *stack = mcEvent->Stack();
51 if(!stack){
52 Printf("%s%d No Stack available",(char*)__FILE__,__LINE__);
53 return;
54 }
55
56 static Int_t iCount = 0;
57 if(iCount>iMaxPrint)return;
58 Int_t nStack = stack->GetNtrack();
59 if(iLast == 0)iLast = nStack;
60 else if(iLast > nStack)iLast = nStack;
61
62
63 Printf("####################################################################");
64 for(Int_t np = iFirst;np<iLast;++np){
65 TParticle *p = stack->Particle(np);
66 Printf("Nr.%d --- Status %d ---- Mother1 %d Mother2 %d Daughter1 %d Daughter2 %d ",
67 np,p->GetStatusCode(),p->GetMother(0),p->GetMother(1),p->GetDaughter(0),p->GetDaughter(1));
68 Printf("Eta %3.3f Phi %3.3f ",p->Eta(),p->Phi());
69 p->Print();
70 Printf("---------------------------------------");
71 }
72 iCount++;
73}
74
75
db6bcb0e 76
77
78void AliAnalysisHelperJetTasks::GetClosestJets(AliAODJet *genJets,const Int_t &nGenJets,
79 AliAODJet *recJets,const Int_t &nRecJets,
80 Int_t *iGenIndex,Int_t *iRecIndex,
81 Int_t iDebug,Float_t maxDist){
82
83 //
84 // Relate the two input jet Arrays
85 //
86
87 //
88 // The association has to be unique
89 // So check in two directions
90 // find the closest rec to a gen
91 // and check if there is no other rec which is closer
92 // Caveat: Close low energy/split jets may disturb this correlation
93
94 // Idea: search in two directions generated e.g (a--e) and rec (1--3)
95 // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
96 // in the end we have something like this
97 // 1 2 3
98 // ------------
99 // a| 3 2 0
100 // b| 0 1 0
101 // c| 0 0 3
102 // d| 0 0 1
103 // e| 0 0 1
104 // Topology
105 // 1 2
106 // a b
107 //
108 // d c
109 // 3 e
110 // Only entries with "3" match from both sides
111
112 const int kMode = 3;
113
114
115
116 if(nRecJets==0||nGenJets==0)return;
117
118 for(int i = 0;i < nGenJets;++i)iGenIndex[i] = -1;
119 for(int j = 0;j < nRecJets;++j)iRecIndex[j] = -1;
120
121 UShort_t *iFlag = new UShort_t[nGenJets*nRecJets];
122 for(int i = 0;i < nGenJets;++i){
123 for(int j = 0;j < nRecJets;++j){
124 iFlag[i*nGenJets+j] = 0;
125 }
126 }
127
128
129
130 // find the closest distance to the generated
131 for(int ig = 0;ig<nGenJets;++ig){
132 Float_t dist = maxDist;
133 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());
134 for(int ir = 0;ir<nRecJets;++ir){
135 Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
136 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());
137 if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
138 if(dR<dist){
139 iRecIndex[ig] = ir;
140 dist = dR;
141 }
142 }
143 if(iRecIndex[ig]>=0)iFlag[ig*nGenJets+iRecIndex[ig]]+=1;
144 // reset...
145 iRecIndex[ig] = -1;
146 }
147 // other way around
148 for(int ir = 0;ir<nRecJets;++ir){
149 Float_t dist = maxDist;
150 for(int ig = 0;ig<nGenJets;++ig){
151 Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
152 if(dR<dist){
153 iGenIndex[ir] = ig;
154 dist = dR;
155 }
156 }
157 if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]*nGenJets+ir]+=2;
158 // reset...
159 iGenIndex[ir] = -1;
160 }
161
162 // check for "true" correlations
163
164 if(iDebug>1)Printf(">>>>>> Matrix");
165
166 for(int ig = 0;ig<nGenJets;++ig){
167 for(int ir = 0;ir<nRecJets;++ir){
168 // Print
169 if(iDebug>1)printf("Flag[%d][%d] %d ",ig,ir,iFlag[ig*nGenJets+ir]);
170
171 if(kMode==3){
172 // we have a uniqie correlation
173 if(iFlag[ig*nGenJets+ir]==3){
174 iGenIndex[ir] = ig;
175 iRecIndex[ig] = ir;
176 }
177 }
178 else{
179 // we just take the correlation from on side
180 if((iFlag[ig*nGenJets+ir]&2)==2){
181 iGenIndex[ir] = ig;
182 }
183 if((iFlag[ig*nGenJets+ir]&1)==1){
184 iRecIndex[ig] = ir;
185 }
186 }
187 }
188 if(iDebug>1)printf("\n");
189 }
190
191 delete [] iFlag;
192
193}
194
195
196
197
198