]>
Commit | Line | Data |
---|---|---|
498224e8 | 1 | //*-- Author: Aleksei Pavlinov (WSU) |
2 | //*-- | |
3 | #include "jetDst.h" | |
4 | #if !defined(__CINT__) || defined(__MAKECINT__) | |
5 | // These files are needed for compilation only | |
6 | // (see directory /auto/u/pavlinov/macros on PDSF) | |
7 | // Using with caution if you are not PAI | |
8 | #include "macroIncludePai.h" | |
9 | #include "macroIncludeAlice.h" | |
30bd7436 | 10 | extern "C++" {void loadlibs();} |
498224e8 | 11 | #endif |
12 | ||
13 | // For convinience only | |
14 | class AliEMCALJetFinder; AliEMCALJetFinder* jetFinder=0; | |
15 | TH2F *hLego, *hLegoW, *hTrig; | |
16 | TH1F *hETrigPatch, *hJetPartPt; | |
17 | TCanvas *c1=0; | |
18 | // directory with file and number of files | |
19 | TString dirIn("/auto/alice/sahal/PPR3/11/"), nFileIn; | |
20 | Int_t file1=1, file2=20, nevMax=0; | |
21 | // change dirOut if you want. | |
22 | TString dirOut("RES/FILE/"), nFileOut; | |
23 | TFile *file=0; | |
24 | // Default value for jet fineer - 11-feb-2202 | |
25 | Int_t debugJet=1; | |
26 | Float_t coneRadius = 0.3; | |
27 | Float_t etSeed = 4.; | |
28 | Float_t minJetEt = 40.; | |
29 | Float_t minCellEt = 1.; // be carefull | |
30 | Float_t minPtCut = 2.; // for including charge particles to file | |
31 | // For BG subtraction | |
32 | Int_t modeSubtraction = 1; | |
33 | Float_t minMove = 0.05; | |
34 | Float_t maxMove = 0.15; | |
35 | Float_t precBg = 0.035; | |
36 | // key 13-feb-2002 | |
37 | Bool_t hadronCorrection = kTRUE; | |
38 | // | |
39 | class AliEMCALJetMicroDst; AliEMCALJetMicroDst *dst=0; | |
40 | Int_t nevSum = 0, nevF; | |
41 | TArrayI nevInFiles(100); | |
42 | ||
498224e8 | 43 | |
44 | void jetDst(Int_t mode, Int_t var, Int_t nmax) | |
45 | { | |
46 | // | |
47 | // mode - variable for selection input mode (see function defineInput()); | |
48 | // var - variable for selection jetfinder mode (see function defineJetFinder()). | |
49 | // | |
50 | // Dynamically link some shared libs | |
51 | if (gClassTable->GetID("AliRun") < 0) { | |
52 | gROOT->LoadMacro("../macros/loadlibs.C"); | |
53 | loadlibs(); | |
54 | } | |
55 | ||
56 | if(!defineInput(mode)) return; | |
57 | ||
58 | dst = new AliEMCALJetMicroDst; | |
59 | dst->SetDebug(1); | |
60 | if(!nFileOut.Contains(".root")){ | |
61 | nFileOut += var; | |
62 | nFileOut += ".root"; | |
63 | } | |
64 | dst->Create(nFileOut.Data()); | |
65 | gROOT->GetListOfBrowsables()->Add(dst); | |
66 | ||
67 | for(Int_t j=file1; j<=file2; j++){ | |
68 | nevF = 0; | |
69 | if(file2>1) {// multiple file | |
70 | nFileIn = dirIn; | |
71 | if (mode>=11 && mode<=11) { | |
72 | nFileIn += j; | |
73 | nFileIn += ("/galice.root"); | |
74 | } else if(mode>=21 && mode<=51){ | |
75 | nFileIn += "galice"; | |
76 | if(j>0) nFileIn += j; | |
77 | nFileIn += ".root"; | |
78 | } | |
79 | } | |
80 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
81 | file = new TFile(nFileIn.Data(), "READ"); | |
82 | printf("\n FILE %i -> %s\n", j, nFileIn.Data()); | |
83 | ||
84 | if (file == 0) continue; | |
85 | ||
86 | // Get AliRun object from file or create it if not on file | |
87 | gAlice = 0; | |
88 | gAlice = (AliRun*)(file->Get("gAlice")); | |
89 | if(!gAlice) return; | |
90 | printf("AliRun object found on file\n"); | |
91 | ||
92 | defineJetFinder(mode,var); | |
93 | // | |
94 | // Loop over events in file | |
95 | // | |
96 | ||
97 | // if(!c1) c1 = new TCanvas("c1","Canvas 1",400,10,600,700); | |
98 | Int_t nhit=0; | |
99 | TString sw; | |
100 | if(nmax>0) nevMax = nmax; | |
101 | for (Int_t nev = 0; nev<= nevMax; nev++) { | |
102 | jetFinder->ResetJets(); | |
103 | Int_t nparticles = gAlice->GetEvent(nev); | |
104 | if (nparticles <= 0) break; // last event | |
105 | nevSum++; | |
106 | nevF++; | |
107 | // Partons jet | |
108 | /* | |
109 | jetFinder->FillFromPartons(); | |
110 | jetFinder->Find(); | |
111 | for(Int_t ij=0; ij<jetFinder->Njets(); ij++){ | |
112 | hJetPartPt->Fill(jetFinder->JetEnergy(ij)); | |
113 | } | |
114 | jetFinder->ResetJets(); | |
115 | */ | |
116 | // ECAL information | |
117 | jetFinder->FillFromHits(); | |
118 | ||
119 | // gROOT->cd(); be carefull with this staff | |
120 | hLegoW = new TH2F((*jetFinder->GetLego())); // for saving | |
121 | sw = hLegoW->GetName(); | |
122 | sw += nev; | |
123 | hLegoW->SetName(sw.Data()); | |
124 | sw = "Energy in ECAL:Event"; | |
125 | sw += nev; | |
126 | hLegoW->SetTitle(sw.Data()); | |
127 | ||
128 | if(hTrig==0) { | |
129 | hLego = jetFinder->GetLego(); | |
130 | hTrig = bookTrigHist(hLego); | |
131 | } | |
132 | fillTriggerPatch(hLegoW, hTrig, hETrigPatch); | |
133 | ||
134 | // jetFinder->FillFromDigits(); | |
135 | // TPC information | |
136 | jetFinder->FillFromTracks(1, 0); | |
137 | // ^ preserves info from hit | |
138 | ||
139 | // TPC information from Hits associated to the EMCAL | |
140 | // jetFinder->FillFromHitFlaggedTracks(1); | |
141 | // | |
142 | jetFinder->Find(); | |
143 | dst->Fill(gAlice, jetFinder); | |
144 | // | |
145 | // Look at results | |
146 | printf("\n Events : %d Jets : %d \n", nev, jetFinder->Njets()); | |
147 | Int_t njet = jetFinder->Njets(); | |
148 | for (Int_t nj=0; nj<njet; nj++) | |
149 | { | |
150 | printf("\n Jet Energy %5d %8.2f \n", | |
151 | nj, jetFinder->JetEnergy(nj)); | |
152 | printf("Jet Phi %5d %8.2f %8.2f \n", | |
153 | nj, jetFinder->JetPhiL(nj), jetFinder->JetPhiW(nj)); | |
154 | printf("Jet Eta %5d %8.2f %8.2f \n", | |
155 | nj, jetFinder->JetEtaL(nj), jetFinder->JetEtaW(nj)); | |
156 | } | |
157 | ||
158 | Int_t in; | |
159 | if (c1 && nev <-10 && njet==1 && jetFinder->JetEnergy(0)>150) { | |
160 | c1->cd(); | |
161 | cout << "Event "<< nev <<endl; | |
162 | hLego->Draw("lego"); | |
163 | c1->Update(); | |
164 | cin >> in; | |
165 | if(in<=0) return; | |
166 | } | |
167 | } // event loop | |
168 | nevInFiles[j-1] = nevF; | |
169 | if (gAlice) { | |
170 | delete gAlice; // don't delete after closing file | |
171 | gAlice = 0; | |
172 | } | |
173 | file->Close(); | |
174 | if(c1) hLego->Draw("lego"); | |
175 | } // file loop | |
176 | if(nmax==0) dst->Close(); | |
177 | nevInFiles.Set(file2); | |
178 | printf("\n Summ.events %6i -> %6i\n", nevSum, (Int_t) nevInFiles.GetSum()); | |
179 | ||
180 | printf("Events "); | |
181 | for(Int_t j=file1; j<=file2; j++){ | |
182 | printf(" %i4 ",nevInFiles[j-1]); | |
183 | if((j-file1)%10==0 && j!=file1) printf("\n "); | |
184 | } | |
185 | } | |
186 | ||
187 | // ################################################################### | |
188 | TString* defineInput(Int_t mode) | |
189 | {// 7-feb-2002 | |
190 | Char_t ctmp[200]; | |
191 | switch(mode){ | |
192 | // case -5 needing for testing after merging | |
193 | case -5: | |
194 | case -4: | |
195 | hadronCorrection = kFALSE; // for checking Tom's idea | |
196 | case -3: | |
197 | case -2: | |
198 | case -1: | |
199 | if (mode==-2) minPtCut = 1.; | |
200 | else if(mode<=-3 || mode>=-5) { | |
201 | minCellEt = 0.0; // 13-feb-2002 | |
202 | minPtCut = 0.0; | |
203 | } | |
204 | ||
205 | coneRadius = 0.5; | |
206 | minJetEt = 20.0; | |
207 | file1=1; | |
208 | file2=1; | |
209 | modeSubtraction = 0; | |
210 | nevMax = 1; | |
211 | nevMax = 999; | |
212 | dirIn = "."; | |
213 | nFileIn = "galice.root"; | |
214 | sprintf(ctmp,"PythiaMode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root", | |
215 | mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut); | |
216 | nFileOut = dirOut + ctmp; | |
217 | if(!hadronCorrection) nFileOut.ReplaceAll(".root","_NoHadrCorr.root"); | |
218 | break; | |
219 | case 1: | |
220 | nFileOut = dirOut + "dstTest.root"; | |
221 | break; | |
222 | ||
223 | case 11: | |
224 | // jj100 41 events per file - see Log file | |
225 | dirIn = "/auto/alice/sahal/PPR4/1/"; | |
226 | file1 = 1; | |
227 | file2 = 30; // max 75 | |
228 | nevMax = 41; | |
229 | sprintf(ctmp,"HijingMode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root", | |
230 | mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut); | |
231 | nFileOut = dirOut + ctmp; | |
232 | break; | |
233 | //============================================================= | |
234 | // central events with jet trigger - 15 feb 2002 | |
235 | //============================================================= | |
236 | case 21: | |
237 | dirIn = "/auto/alice/sahal/kHijing_jj50/"; | |
238 | file1 = 0; | |
239 | file2 = 9; | |
240 | nevMax = 41; | |
241 | sprintf(ctmp,"Hijing_jj50Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root", | |
242 | mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut); | |
243 | nFileOut = dirOut + ctmp; | |
244 | break; | |
245 | case 22: | |
246 | dirIn = "/auto/alice/sahal/kHijing_jj75/"; | |
247 | file1 = 0; | |
248 | file2 = 9; | |
249 | nevMax = 41; | |
250 | sprintf(ctmp,"Hijing_jj75Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root", | |
251 | mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut); | |
252 | nFileOut = dirOut + ctmp; | |
253 | break; | |
254 | case 23: | |
255 | dirIn = "/auto/alice/sahal/kHijing_jj100/"; | |
256 | file1 = 0; | |
257 | file2 = 4; | |
258 | nevMax = 41; | |
259 | sprintf(ctmp,"Hijing_jj100Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root", | |
260 | mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut); | |
261 | nFileOut = dirOut + ctmp; | |
262 | break; | |
263 | case 24: | |
264 | dirIn = "/auto/alice/sahal/kHijing_jj125/"; | |
265 | file1 = 0; | |
266 | file2 = 4; | |
267 | nevMax = 41; | |
268 | sprintf(ctmp,"Hijing_jj125Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root", | |
269 | mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut); | |
270 | nFileOut = dirOut + ctmp; | |
271 | break; | |
272 | ||
273 | case 51: | |
274 | // BG - HIJING events | |
275 | dirIn = "/auto/alice/sahal/background/"; | |
276 | file1 = 1; | |
277 | file2 = 9; | |
278 | nevMax = 100; // ?? | |
279 | sprintf(ctmp,"HijingBGMode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root", | |
280 | mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut); | |
281 | nFileOut = dirOut + ctmp; | |
282 | break; | |
283 | default: | |
284 | printf("\n<I> Wrong input mode => %i\ncout", mode); | |
285 | return 0; | |
286 | } | |
287 | printf("<I> dir %s file1 %3i file2 %3i\n", dirIn.Data(), file1, file2); | |
288 | return &nFileOut; | |
289 | } | |
290 | ||
291 | void defineJetFinder(Int_t mode, Int_t var) | |
292 | {// 11-feb-2002 | |
293 | // Create and configure jet finder | |
294 | // gAlice must be defined already !!! | |
295 | if(!jetFinder) { | |
296 | jetFinder = new AliEMCALJetFinder("UA1 Jet Finder", "Pre Production"); | |
297 | jetFinder->Init(); | |
298 | ||
299 | jetFinder->SetDebug(debugJet); | |
300 | jetFinder->SetConeRadius(coneRadius); | |
301 | jetFinder->SetEtSeed(etSeed); | |
302 | jetFinder->SetMinJetEt(minJetEt); | |
303 | jetFinder->SetMinCellEt(minCellEt); | |
304 | jetFinder->SetPtCut(minPtCut); | |
305 | ||
306 | if(var>=2) jetFinder->SetMomentumSmearing(1); | |
307 | if(var>=3) jetFinder->SetEfficiencySim(1); | |
308 | ||
309 | jetFinder->SetHadronCorrection(0); | |
310 | if(var>=4 && hadronCorrection) { | |
311 | jetFinder->SetHadronCorrection(1); | |
312 | jetFinder->SetHadronCorrector(AliEMCALHadronCorrectionv0::Instance()); | |
313 | } | |
314 | if(var>=5) { // 11-feb-2002 | |
315 | jetFinder->SetParametersForBgSubtraction | |
316 | (modeSubtraction,minMove,maxMove,precBg); | |
317 | } | |
318 | ||
319 | if(mode<0) jetFinder->SetSamplingFraction(1.);//must use only for PYTHIA file | |
320 | else jetFinder->SetSamplingFraction(12.9); | |
321 | gROOT->GetListOfBrowsables()->Add(jetFinder); | |
322 | } | |
323 | } | |
324 | ||
325 | void testDST(char* fname) | |
326 | {// 8-feb-2002 | |
327 | dst = new AliEMCALJetMicroDst; | |
328 | dst->Open(fname); | |
329 | dst->Test(); | |
330 | gROOT->GetListOfBrowsables()->Add(dst); | |
331 | } | |
332 | ||
333 | //----------------- 4-feb-2002 -------------- | |
334 | TH2F* bookTrigHist(TH2F *hid) | |
335 | { | |
336 | if(hid==0) return 0; | |
337 | TAxis *xax = hid->GetXaxis(); | |
338 | Double_t xmin = xax->GetXmin(); | |
339 | Double_t xmax = xax->GetXmax(); | |
340 | Int_t nx = 6; | |
341 | TAxis *yax = hid->GetYaxis(); | |
342 | Double_t ymin = yax->GetXmin(); | |
343 | Double_t ymax = yax->GetXmax(); | |
344 | Int_t ny = 6; | |
345 | TH2F *hout = new TH2F("hTrig","Energy grid for trigger patch", | |
346 | nx,xmin,xmax, ny,ymin,ymax); | |
347 | hout->AddDirectory(0); | |
348 | ||
349 | // time solution | |
350 | hETrigPatch = new TH1F("hETrigPatch","Energy in trigger patch",300, 0.0, 300.); | |
351 | hJetPartPt = new TH1F("hJetPartPt", "E_{T} for partons jet", 300, 0.0, 300.); | |
352 | hETrigPatch->AddDirectory(0); | |
353 | return hout; | |
354 | } | |
355 | ||
356 | void fillTriggerPatch(TH2F *hid, TH2F *htrig, TH1F *hE) | |
357 | { | |
358 | TAxis *xax = hid->GetXaxis(); | |
359 | TAxis *yax = hid->GetYaxis(); | |
360 | Int_t nx = xax->GetNbins(); | |
361 | Int_t ny = yax->GetNbins(); | |
362 | ||
363 | htrig->Reset(); | |
364 | Double_t x, y, e; | |
365 | for(Int_t ix=1; ix<=nx; ix++){ | |
366 | x = xax->GetBinCenter(ix); | |
367 | for(Int_t iy=1; iy<=ny; iy++){ | |
368 | y = yax->GetBinCenter(iy); | |
369 | e = hid->GetBinContent(ix, iy); | |
370 | htrig->Fill(x,y,e); | |
371 | } | |
372 | } | |
373 | ||
374 | xax = htrig->GetXaxis(); | |
375 | yax = htrig->GetYaxis(); | |
376 | nx = xax->GetNbins(); | |
377 | ny = yax->GetNbins(); | |
378 | for(Int_t ix=1; ix<=nx; ix++){ | |
379 | for(Int_t iy=1; iy<=ny; iy++){ | |
380 | hE->Fill(htrig->GetBinContent(ix, iy)); | |
381 | } | |
382 | } | |
383 | } |