]>
Commit | Line | Data |
---|---|---|
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 | ||
15 | ClassImp(AliAnalysisHelperJetTasks) | |
16 | ||
17 | ||
18 | ||
19 | ||
20 | AliGenPythiaEventHeader* 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 | ||
48 | void 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 | ||
78 | void 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 | ||
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 |