]>
Commit | Line | Data |
---|---|---|
f3050824 | 1 | |
2 | #include "TROOT.h" | |
3 | #include "TList.h" | |
188a1ba9 | 4 | #include "TH1F.h" |
5 | #include "TProfile.h" | |
6 | #include "THnSparse.h" | |
7 | #include "TFile.h" | |
f3050824 | 8 | #include "AliMCEvent.h" |
db6bcb0e | 9 | #include "AliAODJet.h" |
f3050824 | 10 | #include "AliStack.h" |
11 | #include "AliGenEventHeader.h" | |
12 | #include "AliGenCocktailEventHeader.h" | |
13 | #include "AliGenPythiaEventHeader.h" | |
14 | #include <fstream> | |
15 | #include <iostream> | |
16 | #include "AliAnalysisHelperJetTasks.h" | |
17 | ||
18 | ||
19 | ClassImp(AliAnalysisHelperJetTasks) | |
20 | ||
21 | ||
22 | ||
23 | ||
24 | AliGenPythiaEventHeader* AliAnalysisHelperJetTasks::GetPythiaEventHeader(AliMCEvent *mcEvent){ | |
25 | ||
26 | AliGenEventHeader* genHeader = mcEvent->GenEventHeader(); | |
27 | AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader); | |
28 | if(!pythiaGenHeader){ | |
29 | // cocktail ?? | |
30 | AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader); | |
31 | ||
32 | if (!genCocktailHeader) { | |
33 | Printf("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__); | |
34 | return 0; | |
35 | } | |
36 | TList* headerList = genCocktailHeader->GetHeaders(); | |
37 | for (Int_t i=0; i<headerList->GetEntries(); i++) { | |
38 | pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i)); | |
39 | if (pythiaGenHeader) | |
40 | break; | |
41 | } | |
42 | if(!pythiaGenHeader){ | |
43 | Printf("%s %d: PythiaHeader not found!",(char*)__FILE__,__LINE__); | |
44 | return 0; | |
45 | } | |
46 | } | |
47 | return pythiaGenHeader; | |
48 | ||
49 | } | |
50 | ||
51 | ||
52 | void AliAnalysisHelperJetTasks::PrintStack(AliMCEvent *mcEvent,Int_t iFirst,Int_t iLast,Int_t iMaxPrint){ | |
53 | ||
54 | AliStack *stack = mcEvent->Stack(); | |
55 | if(!stack){ | |
56 | Printf("%s%d No Stack available",(char*)__FILE__,__LINE__); | |
57 | return; | |
58 | } | |
59 | ||
60 | static Int_t iCount = 0; | |
61 | if(iCount>iMaxPrint)return; | |
62 | Int_t nStack = stack->GetNtrack(); | |
63 | if(iLast == 0)iLast = nStack; | |
64 | else if(iLast > nStack)iLast = nStack; | |
65 | ||
66 | ||
67 | Printf("####################################################################"); | |
68 | for(Int_t np = iFirst;np<iLast;++np){ | |
69 | TParticle *p = stack->Particle(np); | |
70 | Printf("Nr.%d --- Status %d ---- Mother1 %d Mother2 %d Daughter1 %d Daughter2 %d ", | |
71 | np,p->GetStatusCode(),p->GetMother(0),p->GetMother(1),p->GetDaughter(0),p->GetDaughter(1)); | |
72 | Printf("Eta %3.3f Phi %3.3f ",p->Eta(),p->Phi()); | |
73 | p->Print(); | |
74 | Printf("---------------------------------------"); | |
75 | } | |
76 | iCount++; | |
77 | } | |
78 | ||
79 | ||
db6bcb0e | 80 | |
81 | ||
3dc5a1a4 | 82 | void AliAnalysisHelperJetTasks::GetClosestJets(AliAODJet *genJets,const Int_t &kGenJets, |
83 | AliAODJet *recJets,const Int_t &kRecJets, | |
db6bcb0e | 84 | Int_t *iGenIndex,Int_t *iRecIndex, |
85 | Int_t iDebug,Float_t maxDist){ | |
86 | ||
87 | // | |
88 | // Relate the two input jet Arrays | |
89 | // | |
90 | ||
91 | // | |
92 | // The association has to be unique | |
93 | // So check in two directions | |
94 | // find the closest rec to a gen | |
95 | // and check if there is no other rec which is closer | |
96 | // Caveat: Close low energy/split jets may disturb this correlation | |
97 | ||
3dc5a1a4 | 98 | |
db6bcb0e | 99 | // Idea: search in two directions generated e.g (a--e) and rec (1--3) |
100 | // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet | |
101 | // in the end we have something like this | |
102 | // 1 2 3 | |
103 | // ------------ | |
104 | // a| 3 2 0 | |
105 | // b| 0 1 0 | |
106 | // c| 0 0 3 | |
107 | // d| 0 0 1 | |
108 | // e| 0 0 1 | |
109 | // Topology | |
110 | // 1 2 | |
111 | // a b | |
112 | // | |
113 | // d c | |
114 | // 3 e | |
115 | // Only entries with "3" match from both sides | |
3dc5a1a4 | 116 | |
117 | // In case we have more jets than kmaxjets only the | |
118 | // first kmaxjets are searched | |
119 | // all other are -1 | |
120 | // use kMaxJets for a test not to fragemnt the memory... | |
121 | ||
122 | for(int i = 0;i < kGenJets;++i)iGenIndex[i] = -1; | |
123 | for(int j = 0;j < kRecJets;++j)iRecIndex[j] = -1; | |
124 | ||
125 | ||
db6bcb0e | 126 | |
127 | const int kMode = 3; | |
128 | ||
3dc5a1a4 | 129 | const Int_t nGenJets = TMath::Min(kMaxJets,kGenJets); |
130 | const Int_t nRecJets = TMath::Min(kMaxJets,kRecJets); | |
db6bcb0e | 131 | |
132 | if(nRecJets==0||nGenJets==0)return; | |
133 | ||
3dc5a1a4 | 134 | // UShort_t *iFlag = new UShort_t[nGenJets*nRecJets]; |
135 | UShort_t iFlag[kMaxJets*kMaxJets]; | |
db6bcb0e | 136 | for(int i = 0;i < nGenJets;++i){ |
137 | for(int j = 0;j < nRecJets;++j){ | |
138 | iFlag[i*nGenJets+j] = 0; | |
139 | } | |
140 | } | |
141 | ||
142 | ||
143 | ||
144 | // find the closest distance to the generated | |
145 | for(int ig = 0;ig<nGenJets;++ig){ | |
146 | Float_t dist = maxDist; | |
147 | 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()); | |
148 | for(int ir = 0;ir<nRecJets;++ir){ | |
149 | Double_t dR = genJets[ig].DeltaR(&recJets[ir]); | |
150 | 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()); | |
151 | if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR); | |
152 | if(dR<dist){ | |
153 | iRecIndex[ig] = ir; | |
154 | dist = dR; | |
155 | } | |
156 | } | |
157 | if(iRecIndex[ig]>=0)iFlag[ig*nGenJets+iRecIndex[ig]]+=1; | |
158 | // reset... | |
159 | iRecIndex[ig] = -1; | |
160 | } | |
161 | // other way around | |
162 | for(int ir = 0;ir<nRecJets;++ir){ | |
163 | Float_t dist = maxDist; | |
164 | for(int ig = 0;ig<nGenJets;++ig){ | |
165 | Double_t dR = genJets[ig].DeltaR(&recJets[ir]); | |
166 | if(dR<dist){ | |
167 | iGenIndex[ir] = ig; | |
168 | dist = dR; | |
169 | } | |
170 | } | |
171 | if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]*nGenJets+ir]+=2; | |
172 | // reset... | |
173 | iGenIndex[ir] = -1; | |
174 | } | |
175 | ||
176 | // check for "true" correlations | |
177 | ||
178 | if(iDebug>1)Printf(">>>>>> Matrix"); | |
179 | ||
180 | for(int ig = 0;ig<nGenJets;++ig){ | |
181 | for(int ir = 0;ir<nRecJets;++ir){ | |
182 | ||
183 | if(iDebug>1)printf("Flag[%d][%d] %d ",ig,ir,iFlag[ig*nGenJets+ir]); | |
184 | ||
185 | if(kMode==3){ | |
186 | // we have a uniqie correlation | |
187 | if(iFlag[ig*nGenJets+ir]==3){ | |
188 | iGenIndex[ir] = ig; | |
189 | iRecIndex[ig] = ir; | |
190 | } | |
191 | } | |
192 | else{ | |
193 | // we just take the correlation from on side | |
194 | if((iFlag[ig*nGenJets+ir]&2)==2){ | |
195 | iGenIndex[ir] = ig; | |
196 | } | |
197 | if((iFlag[ig*nGenJets+ir]&1)==1){ | |
198 | iRecIndex[ig] = ir; | |
199 | } | |
200 | } | |
201 | } | |
202 | if(iDebug>1)printf("\n"); | |
203 | } | |
db6bcb0e | 204 | } |
205 | ||
206 | ||
207 | ||
188a1ba9 | 208 | void AliAnalysisHelperJetTasks::MergeOutput(char* cFiles, char* cList){ |
db6bcb0e | 209 | |
188a1ba9 | 210 | // This is used to merge the analysis-output from different |
211 | // data samples/pt_hard bins | |
212 | // in case the eventweigth was set to xsection/ntrials already, this | |
213 | // is not needed. Both methods only work in case we do not mix different | |
214 | // pt_hard bins, and do not have overlapping bins | |
db6bcb0e | 215 | |
188a1ba9 | 216 | const Int_t nMaxBins = 12; |
217 | // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05 | |
218 | // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03 | |
219 | // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02 | |
220 | // const Float_t xsection[nBins] = {2.168291,2.44068e-02}; | |
221 | ||
222 | Float_t xsection[nMaxBins]; | |
223 | Float_t nTrials[nMaxBins]; | |
224 | Float_t sf[nMaxBins]; | |
225 | TList *lIn[nMaxBins]; | |
226 | TFile *fIn[nMaxBins]; | |
227 | ||
228 | ifstream in1; | |
229 | in1.open(cFiles); | |
230 | ||
231 | char cFile[120]; | |
232 | Int_t ibTotal = 0; | |
233 | while(in1>>cFile){ | |
234 | fIn[ibTotal] = TFile::Open(cFile); | |
235 | lIn[ibTotal] = (TList*)fIn[ibTotal]->Get(cList); | |
236 | if(!lIn[ibTotal]){ | |
237 | Printf("%s:%d No list %s found, exiting...",__FILE__,__LINE__,cList); | |
238 | fIn[ibTotal]->ls(); | |
239 | return; | |
240 | } | |
241 | TH1* hTrials = (TH1F*)lIn[ibTotal]->FindObject("fh1Trials"); | |
242 | if(!hTrials){ | |
243 | Printf("%s:%d fh1PtHard_Trials not found in list, exiting...",__FILE__,__LINE__); | |
244 | return; | |
245 | } | |
246 | TProfile* hXsec = (TProfile*)lIn[ibTotal]->FindObject("fh1Xsec"); | |
247 | if(!hXsec){ | |
248 | Printf("%s:%d fh1Xsec not found in list, exiting...",__FILE__,__LINE__); | |
249 | return; | |
250 | } | |
251 | xsection[ibTotal] = hXsec->GetBinContent(1); | |
252 | nTrials[ibTotal] = hTrials->Integral(); | |
253 | sf[ibTotal] = xsection[ibTotal]/ nTrials[ibTotal]; | |
254 | ibTotal++; | |
255 | } | |
256 | ||
257 | if(ibTotal==0){ | |
258 | Printf("%s:%d No files found for mergin, exiting",__FILE__,__LINE__); | |
259 | return; | |
260 | } | |
261 | ||
262 | TFile *fOut = new TFile("allpt.root","RECREATE"); | |
263 | TList *lOut = new TList(); | |
264 | lOut->SetName(lIn[0]->GetName()); | |
265 | // for the start scale all... | |
266 | for(int ie = 0; ie < lIn[0]->GetEntries();++ie){ | |
267 | TH1 *h1Add = 0; | |
268 | THnSparse *hnAdd = 0; | |
269 | for(int ib = 0;ib < ibTotal;++ib){ | |
270 | // dynamic cast does not work with cint | |
271 | TObject *h = lIn[ib]->At(ie); | |
272 | if(h->InheritsFrom("TH1")){ | |
273 | TH1 *h1 = (TH1*)h; | |
274 | if(ib==0){ | |
275 | h1Add = (TH1*)h1->Clone(h1->GetName()); | |
276 | h1Add->Scale(sf[ib]); | |
277 | } | |
278 | else{ | |
279 | h1Add->Add(h1,sf[ib]); | |
280 | } | |
281 | } | |
282 | else if(h->InheritsFrom("THnSparse")){ | |
283 | THnSparse *hn = (THnSparse*)h; | |
284 | if(ib==0){ | |
285 | hnAdd = (THnSparse*)hn->Clone(hn->GetName()); | |
286 | hnAdd->Scale(sf[ib]); | |
287 | } | |
288 | else{ | |
289 | hnAdd->Add(hn,sf[ib]); | |
290 | } | |
291 | } | |
292 | ||
293 | ||
294 | }// ib | |
295 | if(h1Add)lOut->Add(h1Add); | |
296 | else if(hnAdd)lOut->Add(hnAdd); | |
297 | } | |
298 | fOut->cd(); | |
299 | lOut->Write(lOut->GetName(),TObject::kSingleKey); | |
300 | fOut->Close(); | |
301 | } |