]>
Commit | Line | Data |
---|---|---|
f3050824 | 1 | |
2 | #include "TROOT.h" | |
519378fb | 3 | #include "TKey.h" |
f3050824 | 4 | #include "TList.h" |
519378fb | 5 | #include "TSystem.h" |
188a1ba9 | 6 | #include "TH1F.h" |
7 | #include "TProfile.h" | |
8 | #include "THnSparse.h" | |
9 | #include "TFile.h" | |
519378fb | 10 | #include "TString.h" |
f3050824 | 11 | #include "AliMCEvent.h" |
5c047edc | 12 | #include "AliLog.h" |
db6bcb0e | 13 | #include "AliAODJet.h" |
f3050824 | 14 | #include "AliStack.h" |
15 | #include "AliGenEventHeader.h" | |
16 | #include "AliGenCocktailEventHeader.h" | |
17 | #include "AliGenPythiaEventHeader.h" | |
18 | #include <fstream> | |
19 | #include <iostream> | |
20 | #include "AliAnalysisHelperJetTasks.h" | |
21 | ||
22 | ||
23 | ClassImp(AliAnalysisHelperJetTasks) | |
24 | ||
25 | ||
26 | ||
27 | ||
28 | AliGenPythiaEventHeader* AliAnalysisHelperJetTasks::GetPythiaEventHeader(AliMCEvent *mcEvent){ | |
29 | ||
30 | AliGenEventHeader* genHeader = mcEvent->GenEventHeader(); | |
31 | AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader); | |
32 | if(!pythiaGenHeader){ | |
33 | // cocktail ?? | |
34 | AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader); | |
35 | ||
36 | if (!genCocktailHeader) { | |
5c047edc | 37 | AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)"); |
38 | // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__)); | |
f3050824 | 39 | return 0; |
40 | } | |
41 | TList* headerList = genCocktailHeader->GetHeaders(); | |
42 | for (Int_t i=0; i<headerList->GetEntries(); i++) { | |
43 | pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i)); | |
44 | if (pythiaGenHeader) | |
45 | break; | |
46 | } | |
47 | if(!pythiaGenHeader){ | |
5c047edc | 48 | AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found"); |
f3050824 | 49 | return 0; |
50 | } | |
51 | } | |
52 | return pythiaGenHeader; | |
53 | ||
54 | } | |
55 | ||
56 | ||
57 | void AliAnalysisHelperJetTasks::PrintStack(AliMCEvent *mcEvent,Int_t iFirst,Int_t iLast,Int_t iMaxPrint){ | |
58 | ||
59 | AliStack *stack = mcEvent->Stack(); | |
60 | if(!stack){ | |
61 | Printf("%s%d No Stack available",(char*)__FILE__,__LINE__); | |
62 | return; | |
63 | } | |
64 | ||
65 | static Int_t iCount = 0; | |
66 | if(iCount>iMaxPrint)return; | |
67 | Int_t nStack = stack->GetNtrack(); | |
68 | if(iLast == 0)iLast = nStack; | |
69 | else if(iLast > nStack)iLast = nStack; | |
70 | ||
71 | ||
72 | Printf("####################################################################"); | |
73 | for(Int_t np = iFirst;np<iLast;++np){ | |
74 | TParticle *p = stack->Particle(np); | |
75 | Printf("Nr.%d --- Status %d ---- Mother1 %d Mother2 %d Daughter1 %d Daughter2 %d ", | |
76 | np,p->GetStatusCode(),p->GetMother(0),p->GetMother(1),p->GetDaughter(0),p->GetDaughter(1)); | |
77 | Printf("Eta %3.3f Phi %3.3f ",p->Eta(),p->Phi()); | |
78 | p->Print(); | |
79 | Printf("---------------------------------------"); | |
80 | } | |
81 | iCount++; | |
82 | } | |
83 | ||
84 | ||
db6bcb0e | 85 | |
86 | ||
3dc5a1a4 | 87 | void AliAnalysisHelperJetTasks::GetClosestJets(AliAODJet *genJets,const Int_t &kGenJets, |
88 | AliAODJet *recJets,const Int_t &kRecJets, | |
db6bcb0e | 89 | Int_t *iGenIndex,Int_t *iRecIndex, |
90 | Int_t iDebug,Float_t maxDist){ | |
91 | ||
92 | // | |
93 | // Relate the two input jet Arrays | |
94 | // | |
95 | ||
96 | // | |
97 | // The association has to be unique | |
98 | // So check in two directions | |
99 | // find the closest rec to a gen | |
100 | // and check if there is no other rec which is closer | |
101 | // Caveat: Close low energy/split jets may disturb this correlation | |
102 | ||
3dc5a1a4 | 103 | |
db6bcb0e | 104 | // Idea: search in two directions generated e.g (a--e) and rec (1--3) |
105 | // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet | |
106 | // in the end we have something like this | |
107 | // 1 2 3 | |
108 | // ------------ | |
109 | // a| 3 2 0 | |
110 | // b| 0 1 0 | |
111 | // c| 0 0 3 | |
112 | // d| 0 0 1 | |
113 | // e| 0 0 1 | |
114 | // Topology | |
115 | // 1 2 | |
116 | // a b | |
117 | // | |
118 | // d c | |
119 | // 3 e | |
120 | // Only entries with "3" match from both sides | |
3dc5a1a4 | 121 | |
122 | // In case we have more jets than kmaxjets only the | |
123 | // first kmaxjets are searched | |
124 | // all other are -1 | |
125 | // use kMaxJets for a test not to fragemnt the memory... | |
126 | ||
127 | for(int i = 0;i < kGenJets;++i)iGenIndex[i] = -1; | |
128 | for(int j = 0;j < kRecJets;++j)iRecIndex[j] = -1; | |
129 | ||
130 | ||
db6bcb0e | 131 | |
132 | const int kMode = 3; | |
133 | ||
3dc5a1a4 | 134 | const Int_t nGenJets = TMath::Min(kMaxJets,kGenJets); |
135 | const Int_t nRecJets = TMath::Min(kMaxJets,kRecJets); | |
db6bcb0e | 136 | |
137 | if(nRecJets==0||nGenJets==0)return; | |
138 | ||
3dc5a1a4 | 139 | // UShort_t *iFlag = new UShort_t[nGenJets*nRecJets]; |
140 | UShort_t iFlag[kMaxJets*kMaxJets]; | |
db6bcb0e | 141 | for(int i = 0;i < nGenJets;++i){ |
142 | for(int j = 0;j < nRecJets;++j){ | |
143 | iFlag[i*nGenJets+j] = 0; | |
144 | } | |
145 | } | |
146 | ||
147 | ||
148 | ||
149 | // find the closest distance to the generated | |
150 | for(int ig = 0;ig<nGenJets;++ig){ | |
151 | Float_t dist = maxDist; | |
152 | 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()); | |
153 | for(int ir = 0;ir<nRecJets;++ir){ | |
154 | Double_t dR = genJets[ig].DeltaR(&recJets[ir]); | |
155 | 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()); | |
156 | if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR); | |
157 | if(dR<dist){ | |
158 | iRecIndex[ig] = ir; | |
159 | dist = dR; | |
160 | } | |
161 | } | |
162 | if(iRecIndex[ig]>=0)iFlag[ig*nGenJets+iRecIndex[ig]]+=1; | |
163 | // reset... | |
164 | iRecIndex[ig] = -1; | |
165 | } | |
166 | // other way around | |
167 | for(int ir = 0;ir<nRecJets;++ir){ | |
168 | Float_t dist = maxDist; | |
169 | for(int ig = 0;ig<nGenJets;++ig){ | |
170 | Double_t dR = genJets[ig].DeltaR(&recJets[ir]); | |
171 | if(dR<dist){ | |
172 | iGenIndex[ir] = ig; | |
173 | dist = dR; | |
174 | } | |
175 | } | |
176 | if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]*nGenJets+ir]+=2; | |
177 | // reset... | |
178 | iGenIndex[ir] = -1; | |
179 | } | |
180 | ||
181 | // check for "true" correlations | |
182 | ||
183 | if(iDebug>1)Printf(">>>>>> Matrix"); | |
184 | ||
185 | for(int ig = 0;ig<nGenJets;++ig){ | |
186 | for(int ir = 0;ir<nRecJets;++ir){ | |
187 | ||
188 | if(iDebug>1)printf("Flag[%d][%d] %d ",ig,ir,iFlag[ig*nGenJets+ir]); | |
189 | ||
190 | if(kMode==3){ | |
191 | // we have a uniqie correlation | |
192 | if(iFlag[ig*nGenJets+ir]==3){ | |
193 | iGenIndex[ir] = ig; | |
194 | iRecIndex[ig] = ir; | |
195 | } | |
196 | } | |
197 | else{ | |
198 | // we just take the correlation from on side | |
199 | if((iFlag[ig*nGenJets+ir]&2)==2){ | |
200 | iGenIndex[ir] = ig; | |
201 | } | |
202 | if((iFlag[ig*nGenJets+ir]&1)==1){ | |
203 | iRecIndex[ig] = ir; | |
204 | } | |
205 | } | |
206 | } | |
207 | if(iDebug>1)printf("\n"); | |
208 | } | |
db6bcb0e | 209 | } |
210 | ||
211 | ||
212 | ||
188a1ba9 | 213 | void AliAnalysisHelperJetTasks::MergeOutput(char* cFiles, char* cList){ |
db6bcb0e | 214 | |
188a1ba9 | 215 | // This is used to merge the analysis-output from different |
216 | // data samples/pt_hard bins | |
217 | // in case the eventweigth was set to xsection/ntrials already, this | |
218 | // is not needed. Both methods only work in case we do not mix different | |
219 | // pt_hard bins, and do not have overlapping bins | |
db6bcb0e | 220 | |
188a1ba9 | 221 | const Int_t nMaxBins = 12; |
222 | // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05 | |
223 | // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03 | |
224 | // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02 | |
225 | // const Float_t xsection[nBins] = {2.168291,2.44068e-02}; | |
226 | ||
227 | Float_t xsection[nMaxBins]; | |
228 | Float_t nTrials[nMaxBins]; | |
229 | Float_t sf[nMaxBins]; | |
230 | TList *lIn[nMaxBins]; | |
231 | TFile *fIn[nMaxBins]; | |
232 | ||
233 | ifstream in1; | |
234 | in1.open(cFiles); | |
235 | ||
236 | char cFile[120]; | |
237 | Int_t ibTotal = 0; | |
238 | while(in1>>cFile){ | |
239 | fIn[ibTotal] = TFile::Open(cFile); | |
240 | lIn[ibTotal] = (TList*)fIn[ibTotal]->Get(cList); | |
5c047edc | 241 | Printf("Merging file %s",cFile); |
188a1ba9 | 242 | if(!lIn[ibTotal]){ |
243 | Printf("%s:%d No list %s found, exiting...",__FILE__,__LINE__,cList); | |
244 | fIn[ibTotal]->ls(); | |
245 | return; | |
246 | } | |
247 | TH1* hTrials = (TH1F*)lIn[ibTotal]->FindObject("fh1Trials"); | |
248 | if(!hTrials){ | |
249 | Printf("%s:%d fh1PtHard_Trials not found in list, exiting...",__FILE__,__LINE__); | |
250 | return; | |
251 | } | |
252 | TProfile* hXsec = (TProfile*)lIn[ibTotal]->FindObject("fh1Xsec"); | |
253 | if(!hXsec){ | |
254 | Printf("%s:%d fh1Xsec not found in list, exiting...",__FILE__,__LINE__); | |
255 | return; | |
256 | } | |
257 | xsection[ibTotal] = hXsec->GetBinContent(1); | |
258 | nTrials[ibTotal] = hTrials->Integral(); | |
259 | sf[ibTotal] = xsection[ibTotal]/ nTrials[ibTotal]; | |
260 | ibTotal++; | |
261 | } | |
262 | ||
263 | if(ibTotal==0){ | |
264 | Printf("%s:%d No files found for mergin, exiting",__FILE__,__LINE__); | |
265 | return; | |
266 | } | |
267 | ||
268 | TFile *fOut = new TFile("allpt.root","RECREATE"); | |
269 | TList *lOut = new TList(); | |
270 | lOut->SetName(lIn[0]->GetName()); | |
271 | // for the start scale all... | |
272 | for(int ie = 0; ie < lIn[0]->GetEntries();++ie){ | |
273 | TH1 *h1Add = 0; | |
274 | THnSparse *hnAdd = 0; | |
275 | for(int ib = 0;ib < ibTotal;++ib){ | |
276 | // dynamic cast does not work with cint | |
277 | TObject *h = lIn[ib]->At(ie); | |
278 | if(h->InheritsFrom("TH1")){ | |
279 | TH1 *h1 = (TH1*)h; | |
280 | if(ib==0){ | |
281 | h1Add = (TH1*)h1->Clone(h1->GetName()); | |
282 | h1Add->Scale(sf[ib]); | |
283 | } | |
284 | else{ | |
285 | h1Add->Add(h1,sf[ib]); | |
286 | } | |
287 | } | |
288 | else if(h->InheritsFrom("THnSparse")){ | |
289 | THnSparse *hn = (THnSparse*)h; | |
290 | if(ib==0){ | |
291 | hnAdd = (THnSparse*)hn->Clone(hn->GetName()); | |
292 | hnAdd->Scale(sf[ib]); | |
293 | } | |
294 | else{ | |
295 | hnAdd->Add(hn,sf[ib]); | |
296 | } | |
297 | } | |
298 | ||
299 | ||
300 | }// ib | |
301 | if(h1Add)lOut->Add(h1Add); | |
302 | else if(hnAdd)lOut->Add(hnAdd); | |
303 | } | |
304 | fOut->cd(); | |
305 | lOut->Write(lOut->GetName(),TObject::kSingleKey); | |
306 | fOut->Close(); | |
307 | } | |
519378fb | 308 | |
309 | Bool_t AliAnalysisHelperJetTasks::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){ | |
310 | // | |
311 | // get the cross section and the trails either from pyxsec.root or from pysec_hists.root | |
312 | // This is to called in Notify and should provide the path to the AOD/ESD file | |
313 | ||
314 | TString file(currFile); | |
315 | fXsec = 0; | |
316 | fTrials = 1; | |
317 | ||
318 | if(file.Contains("root_archive.zip#")){ | |
319 | Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact); | |
320 | Ssiz_t pos = file.Index("#",1,pos1,TString::kExact); | |
321 | file.Replace(pos+1,20,""); | |
322 | } | |
323 | else { | |
324 | // not an archive take the basename.... | |
325 | file.ReplaceAll(gSystem->BaseName(file.Data()),""); | |
326 | } | |
327 | Printf("%s",file.Data()); | |
328 | ||
329 | ||
330 | ||
331 | ||
332 | TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root | |
333 | if(!fxsec){ | |
334 | // next trial fetch the histgram file | |
335 | fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root")); | |
336 | if(!fxsec){ | |
337 | // not a severe condition but inciate that we have no information | |
338 | return kFALSE; | |
339 | } | |
340 | else{ | |
341 | // find the tlist we want to be independtent of the name so use the Tkey | |
342 | TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); | |
343 | if(!key){ | |
344 | fxsec->Close(); | |
345 | return kFALSE; | |
346 | } | |
347 | TList *list = dynamic_cast<TList*>(key->ReadObj()); | |
348 | if(!list){ | |
349 | fxsec->Close(); | |
350 | return kFALSE; | |
351 | } | |
352 | fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1); | |
353 | fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1); | |
354 | fxsec->Close(); | |
355 | } | |
356 | } // no tree pyxsec.root | |
357 | else { | |
358 | TTree *xtree = (TTree*)fxsec->Get("Xsection"); | |
359 | if(!xtree){ | |
360 | fxsec->Close(); | |
361 | return kFALSE; | |
362 | } | |
363 | UInt_t ntrials = 0; | |
364 | Double_t xsection = 0; | |
365 | xtree->SetBranchAddress("xsection",&xsection); | |
366 | xtree->SetBranchAddress("ntrials",&ntrials); | |
367 | xtree->GetEntry(0); | |
368 | fTrials = ntrials; | |
369 | fXsec = xsection; | |
370 | fxsec->Close(); | |
371 | } | |
372 | return kTRUE; | |
373 | } |