]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EMCAL/jetDst.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Debug printout commented out
[u/mrichter/AliRoot.git] / EMCAL / jetDst.C
CommitLineData
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 10extern "C++" {void loadlibs();}
498224e8 11#endif
12
13// For convinience only
14class AliEMCALJetFinder; AliEMCALJetFinder* jetFinder=0;
15TH2F *hLego, *hLegoW, *hTrig;
16TH1F *hETrigPatch, *hJetPartPt;
17TCanvas *c1=0;
18// directory with file and number of files
19TString dirIn("/auto/alice/sahal/PPR3/11/"), nFileIn;
20Int_t file1=1, file2=20, nevMax=0;
21// change dirOut if you want.
22TString dirOut("RES/FILE/"), nFileOut;
23TFile *file=0;
24// Default value for jet fineer - 11-feb-2202
25Int_t debugJet=1;
26Float_t coneRadius = 0.3;
27Float_t etSeed = 4.;
28Float_t minJetEt = 40.;
29Float_t minCellEt = 1.; // be carefull
30Float_t minPtCut = 2.; // for including charge particles to file
31// For BG subtraction
32Int_t modeSubtraction = 1;
33Float_t minMove = 0.05;
34Float_t maxMove = 0.15;
35Float_t precBg = 0.035;
36// key 13-feb-2002
37Bool_t hadronCorrection = kTRUE;
38//
39class AliEMCALJetMicroDst; AliEMCALJetMicroDst *dst=0;
40Int_t nevSum = 0, nevF;
41TArrayI nevInFiles(100);
42
498224e8 43
44void 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// ###################################################################
188TString* 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
291void 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
325void 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 --------------
334TH2F* 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
356void 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}