1 |
mocchiut |
1.1 |
#include <stdlib.h> |
2 |
|
|
#include <iostream> |
3 |
|
|
|
4 |
|
|
#include <TString.h> |
5 |
|
|
#include <TH1F.h> |
6 |
|
|
#include <TH2F.h> |
7 |
|
|
#include <TMath.h> |
8 |
|
|
#include <TLine.h> |
9 |
|
|
#include <TPolyMarker.h> |
10 |
|
|
#include <TSelector.h> |
11 |
|
|
#include <TFile.h> |
12 |
|
|
#include <TGraphErrors.h> |
13 |
|
|
#include <TMultiGraph.h> |
14 |
|
|
#include <TCanvas.h> |
15 |
|
|
|
16 |
|
|
#include <PamLevel2.h> |
17 |
|
|
#include <CaloCore.h> |
18 |
|
|
#include <CaloREPCore.h> |
19 |
|
|
//#include <caloclassesfun.h> |
20 |
|
|
|
21 |
|
|
using namespace std; |
22 |
|
|
|
23 |
|
|
void CaloREPCore(TString ddir,TString list, TString options, TString outfile){ |
24 |
|
|
// |
25 |
|
|
// |
26 |
|
|
// |
27 |
|
|
Bool_t verbose = true; |
28 |
|
|
Bool_t debug = false; |
29 |
|
|
// |
30 |
|
|
Bool_t trackanyway = true; |
31 |
|
|
// |
32 |
|
|
Float_t rigdefault = 50.; |
33 |
|
|
// |
34 |
|
|
Bool_t hZn = true; |
35 |
|
|
// |
36 |
|
|
Bool_t withtrk = true; |
37 |
|
|
// |
38 |
|
|
Bool_t st = true; |
39 |
|
|
// |
40 |
|
|
Bool_t mechal = false; |
41 |
|
|
// |
42 |
|
|
Float_t tmptrigty = -1.; |
43 |
|
|
Int_t S3 = 0; |
44 |
|
|
Int_t S2 = 0; |
45 |
|
|
Int_t S12 = 0; |
46 |
|
|
Int_t S11 = 0; |
47 |
|
|
TString host = "mysql://localhost/pamelaprod"; |
48 |
|
|
TString user = "anonymous"; |
49 |
|
|
TString psw = ""; |
50 |
|
|
const char *pamdbhost=gSystem->Getenv("PAM_DBHOST"); |
51 |
|
|
const char *pamdbuser=gSystem->Getenv("PAM_DBUSER"); |
52 |
|
|
const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW"); |
53 |
|
|
if ( !pamdbhost ) pamdbhost = ""; |
54 |
|
|
if ( !pamdbuser ) pamdbuser = ""; |
55 |
|
|
if ( !pamdbpsw ) pamdbpsw = ""; |
56 |
|
|
if ( strcmp(pamdbhost,"") ) host = pamdbhost; |
57 |
|
|
if ( strcmp(pamdbuser,"") ) user = pamdbuser; |
58 |
|
|
if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw; |
59 |
|
|
GL_TABLES *glt = new GL_TABLES(host,user,psw); |
60 |
|
|
// -------------------- |
61 |
|
|
// read input file/list |
62 |
|
|
// -------------------- |
63 |
|
|
TString outfile_t =0; |
64 |
|
|
if( outfile.IsNull() ){ |
65 |
|
|
if(!list.IsNull())outfile = list(0,list.Last('.')); |
66 |
|
|
else outfile = "output"; |
67 |
|
|
}else{ |
68 |
|
|
if(outfile.Contains(".root"))outfile = outfile(0,outfile.Last('.')); |
69 |
|
|
}; |
70 |
|
|
TFile *outft = NULL; |
71 |
|
|
outfile_t = outfile; |
72 |
|
|
outfile_t.Append("-tree.root"); |
73 |
|
|
outft = (TFile*)gROOT->FindObject(outfile_t); if (outft) outft->Close(); |
74 |
|
|
outft = new TFile(outfile_t,"RECREATE"); |
75 |
|
|
if(outft->IsZombie()){ |
76 |
|
|
cout << "Output file could not be created\n"; |
77 |
|
|
return; |
78 |
|
|
}; |
79 |
|
|
cout << "Created output file: "<<outft->GetName()<<endl; |
80 |
|
|
// |
81 |
|
|
PamLevel2 *event = new PamLevel2(ddir,list,options); |
82 |
|
|
event->SetMaxShift(500); |
83 |
|
|
// |
84 |
|
|
TrkParams::Load(); |
85 |
|
|
// |
86 |
|
|
GL_ROOT *glroot = new GL_ROOT(); |
87 |
|
|
// |
88 |
|
|
CaloLevel1 *c1 = new CaloLevel1(); |
89 |
|
|
CaloLevel2 *ca = new CaloLevel2(); |
90 |
|
|
CaloLevel0 *c0 = new CaloLevel0(); |
91 |
|
|
// |
92 |
|
|
outft->cd(); |
93 |
|
|
event->CreateCloneTrees(outft); |
94 |
|
|
outft->cd(); |
95 |
|
|
TTree *calo = new TTree("Calorimeter","PAMELA Level2 calorimeter data"); |
96 |
|
|
calo->SetAutoSave(900000000000000LL); |
97 |
|
|
ca->Set(); |
98 |
|
|
calo->Branch("CaloLevel2","CaloLevel2",&ca); |
99 |
|
|
calo->Branch("CaloLevel1","CaloLevel1",&c1); |
100 |
|
|
// |
101 |
|
|
// ----------------- |
102 |
|
|
// loop over events |
103 |
|
|
// ----------------- |
104 |
|
|
// |
105 |
|
|
UInt_t srun = 0; |
106 |
|
|
ULong64_t nevents = event->GetEntries(); |
107 |
|
|
// |
108 |
|
|
cout << endl<<" Start loop over events: "<< nevents<<endl; |
109 |
|
|
// |
110 |
|
|
for(ULong64_t iev=0; iev<nevents; iev++){ |
111 |
|
|
//for(ULong64_t iev=0; iev<10000; iev++){ |
112 |
|
|
// |
113 |
|
|
if ( !(iev%1000) ) printf(" %lluK \n",iev/1000); |
114 |
|
|
event->Clear(); |
115 |
|
|
// |
116 |
|
|
if( event->GetEntry(iev) ){ |
117 |
|
|
// |
118 |
|
|
if ( event->GetYodaEntry() ){ |
119 |
|
|
// |
120 |
|
|
event->FillCloneTrees(); |
121 |
|
|
// |
122 |
|
|
Int_t ye = event->GetYodaTree()->GetReadEntry(); |
123 |
|
|
// |
124 |
|
|
UInt_t evfrom = ye; |
125 |
|
|
UInt_t totevent = 1; |
126 |
|
|
Int_t sgnl = 0; |
127 |
|
|
// |
128 |
|
|
// for every new run |
129 |
|
|
// |
130 |
|
|
if ( srun != event->GetRunInfo()->ID ){ |
131 |
|
|
// |
132 |
|
|
UInt_t runheadtime = event->GetRunInfo()->RUNHEADER_TIME; |
133 |
|
|
UInt_t runtrailtime = event->GetRunInfo()->RUNTRAILER_TIME; |
134 |
|
|
srun = event->GetRunInfo()->ID; |
135 |
|
|
// |
136 |
|
|
// prepare the timesync for the db |
137 |
|
|
// |
138 |
|
|
TString host = glt->CGetHost(); |
139 |
|
|
TString user = glt->CGetUser(); |
140 |
|
|
TString psw = glt->CGetPsw(); |
141 |
|
|
TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data()); |
142 |
|
|
if ( !dbc->IsConnected() ) return; |
143 |
|
|
stringstream myquery; |
144 |
|
|
myquery.str(""); |
145 |
|
|
myquery << "SET time_zone='+0:00'"; |
146 |
|
|
dbc->Query(myquery.str().c_str()); |
147 |
|
|
// |
148 |
|
|
// Search in the DB the path and name of the LEVEL0 file to be processed. |
149 |
|
|
// |
150 |
|
|
glroot->Query_GL_ROOT(event->GetRunInfo()->ID_ROOT_L0,dbc); |
151 |
|
|
// |
152 |
|
|
TString fname; |
153 |
|
|
stringstream ftmpname; |
154 |
|
|
ftmpname.str(""); |
155 |
|
|
ftmpname << glroot->PATH.Data() << "/"; |
156 |
|
|
ftmpname << glroot->NAME.Data(); |
157 |
|
|
fname = ftmpname.str().c_str(); |
158 |
|
|
// |
159 |
|
|
// print out informations |
160 |
|
|
// |
161 |
|
|
if ( verbose ) printf("\n LEVEL0 data file: %s \n",fname.Data()); |
162 |
|
|
if ( verbose ) printf(" RUN HEADER absolute time is: %u \n",runheadtime); |
163 |
|
|
if ( verbose ) printf(" RUN TRAILER absolute time is: %u \n",runtrailtime); |
164 |
|
|
if ( verbose ) printf(" Run %u: from %i to %i (runheader %u runtrailer %u )\n\n",event->GetRunInfo()->ID,evfrom,evfrom+totevent,event->GetRunInfo()->EV_FROM,event->GetRunInfo()->EV_TO); |
165 |
|
|
// |
166 |
|
|
// Construct the event object, look for the calibration which include the first header |
167 |
|
|
// |
168 |
|
|
if ( verbose ) printf(" Check for calorimeter calibrations and initialize event object \n"); |
169 |
|
|
// |
170 |
|
|
c0->ProcessingInit(glt,runheadtime,sgnl,event->GetYodaTree(),debug,verbose); |
171 |
|
|
c0->SetCrossTalk(true); |
172 |
|
|
c0->SetCrossTalkType(false); |
173 |
|
|
// c0->SetCrossTalkType(true); |
174 |
|
|
// |
175 |
|
|
// Check if we have to load parameter files |
176 |
|
|
// |
177 |
|
|
sgnl = 0; |
178 |
|
|
sgnl = c0->ChkParam(glt,runheadtime,mechal); // calorimeter parameter files |
179 |
|
|
// |
180 |
|
|
// Calculate the cross talk corrections needed for this run using calibration informations |
181 |
|
|
// |
182 |
|
|
sgnl = 0; |
183 |
|
|
sgnl = c0->CalcCrossTalkCorr(glt,runheadtime); |
184 |
|
|
printf(" sgnl from crosstalkcorr %i \n",sgnl); |
185 |
|
|
// |
186 |
|
|
if ( dbc ){ |
187 |
|
|
dbc->Close(); |
188 |
|
|
delete dbc; |
189 |
|
|
}; |
190 |
|
|
// |
191 |
|
|
}; |
192 |
|
|
// |
193 |
|
|
// printf(" ye is %i \n",ye); |
194 |
|
|
// |
195 |
|
|
// run over all the events of the run |
196 |
|
|
// |
197 |
|
|
UInt_t atime = event->GetOrbitalInfo()->absTime; |
198 |
|
|
// |
199 |
|
|
// start processing |
200 |
|
|
// |
201 |
|
|
c1->Clear(); |
202 |
|
|
ca->Clear(); |
203 |
|
|
// |
204 |
|
|
// determine who generate the trigger for this event (TOF, S4/PULSER, CALO) |
205 |
|
|
// |
206 |
|
|
S3 = 0; |
207 |
|
|
S2 = 0; |
208 |
|
|
S12 = 0; |
209 |
|
|
S11 = 0; |
210 |
|
|
S3 = event->GetTrigLevel2()->patterntrig[2]; |
211 |
|
|
S2 = event->GetTrigLevel2()->patterntrig[3]; |
212 |
|
|
S12 = event->GetTrigLevel2()->patterntrig[4]; |
213 |
|
|
S11 = event->GetTrigLevel2()->patterntrig[5]; |
214 |
|
|
if ( event->GetTrigLevel2()->patterntrig[1] & (1<<0) ) tmptrigty = 1.; |
215 |
|
|
if ( event->GetTrigLevel2()->patterntrig[0] ) tmptrigty = 2.; |
216 |
|
|
if ( S3 || S2 || S12 || S11 ) tmptrigty = 0.; |
217 |
|
|
if ( !(event->GetTrigLevel2()->patterntrig[1] & (1<<0)) && !event->GetTrigLevel2()->patterntrig[0] && !S3 && !S2 && !S12 && !S11 ) tmptrigty = 1.; |
218 |
|
|
c0->clevel2->trigty = tmptrigty; |
219 |
|
|
// |
220 |
|
|
// check if the calibration we are using is still good, if not load another calibration |
221 |
|
|
// |
222 |
|
|
sgnl = 0; |
223 |
|
|
sgnl = c0->ChkCalib(glt,atime); |
224 |
|
|
// |
225 |
|
|
// do we have at least one track from the tracker? this check has been disabled |
226 |
|
|
// |
227 |
|
|
c0->clevel1->good2 = 1; |
228 |
|
|
// |
229 |
|
|
// use the whole calorimeter, 22 W planes |
230 |
|
|
// |
231 |
|
|
c0->clevel1->npla = 22; |
232 |
|
|
// |
233 |
|
|
// Use the standard silicon planes shifting |
234 |
|
|
// |
235 |
|
|
c0->clevel1->reverse = 0; |
236 |
|
|
// |
237 |
|
|
// |
238 |
|
|
// Calibrate calorimeter event "re" and store output in the two structures that will be passed to fortran routine |
239 |
|
|
// |
240 |
|
|
c0->Calibrate(ye); |
241 |
|
|
// |
242 |
|
|
// Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract topological variables. |
243 |
|
|
// |
244 |
|
|
// |
245 |
|
|
// Calculate variables common to all tracks (qtot, nstrip, etc.) |
246 |
|
|
// |
247 |
|
|
c0->GetCommonVar(); |
248 |
|
|
// |
249 |
|
|
// Fill common variables |
250 |
|
|
// |
251 |
|
|
c0->FillCommonVar(c1,ca); |
252 |
|
|
// |
253 |
|
|
// Calculate variables related to tracks only if we have at least one track (from selftrigger and/or tracker) |
254 |
|
|
// |
255 |
|
|
UInt_t ntrkentry = 0; |
256 |
|
|
// |
257 |
|
|
Bool_t filled = false; |
258 |
|
|
// |
259 |
|
|
// Run over tracks (tracker or calorimeter ) |
260 |
|
|
// |
261 |
|
|
if ( withtrk ){ |
262 |
|
|
for(Int_t nt=0; nt < event->GetTrkLevel2()->ntrk(); nt++){ |
263 |
|
|
// |
264 |
|
|
c0->clevel1->good2 = 1; |
265 |
|
|
// |
266 |
|
|
TrkTrack *ptt = event->GetTrkLevel2()->GetStoredTrack(nt)->GetTrkTrack(); |
267 |
|
|
// |
268 |
|
|
c0->clevel1->trkchi2 = 0; |
269 |
|
|
// |
270 |
|
|
// Copy the alpha vector in the input structure |
271 |
|
|
// |
272 |
|
|
for (Int_t e = 0; e < 5 ; e++){ |
273 |
|
|
c0->clevel1->al_p[e][0] = ptt->al[e]; |
274 |
|
|
}; |
275 |
|
|
// |
276 |
|
|
// Get tracker related variables for this track |
277 |
|
|
// |
278 |
|
|
c0->GetTrkVar(); |
279 |
|
|
// |
280 |
|
|
// Save tracker track sequence number |
281 |
|
|
// |
282 |
|
|
c0->trkseqno = nt; |
283 |
|
|
// |
284 |
|
|
// Copy values in the class ca from the structure clevel2 |
285 |
|
|
// |
286 |
|
|
c0->FillTrkVar(ca,ntrkentry); |
287 |
|
|
ntrkentry++; |
288 |
|
|
filled = true; |
289 |
|
|
// |
290 |
|
|
}; // loop on all the tracks |
291 |
|
|
}; |
292 |
|
|
// |
293 |
|
|
// if no tracks found but there is the possibility to have a good track we should try to calculate anyway the track related variables using the calorimeter |
294 |
|
|
// fit of the track (to be used for example when TRK is off due to any reason like IPM3/5 off). |
295 |
|
|
// here we make an event selection so it must be done very carefully... |
296 |
|
|
// |
297 |
|
|
// conditions are: 0) no track from the tracker 1) we have a track fit both in x and y 2) no problems with calo for this event 3) no selftrigger event |
298 |
|
|
// |
299 |
|
|
if ( trackanyway && !filled && c0->clevel2->npcfit[0] >= 2 && c0->clevel2->npcfit[1] >= 2 && c0->clevel2->good != 0 && c0->clevel2->trigty < 2. ){ |
300 |
|
|
if ( debug ) printf(" Event with a track not fitted by the tracker at entry %i \n",ye); |
301 |
|
|
// |
302 |
|
|
// Disable "track mode" in the fortran routine |
303 |
|
|
// |
304 |
|
|
c0->clevel1->good2 = 0; |
305 |
|
|
c0->clevel1->riginput = rigdefault; |
306 |
|
|
// |
307 |
|
|
// We have a selftrigger event to analyze. |
308 |
|
|
// |
309 |
|
|
for (Int_t e = 0; e < 5 ; e++){ |
310 |
|
|
c0->clevel1->al_p[e][0] = 0.; |
311 |
|
|
c0->clevel1->al_p[e][1] = 0.; |
312 |
|
|
}; |
313 |
|
|
c0->clevel1->trkchi2 = 0; |
314 |
|
|
// |
315 |
|
|
c0->GetTrkVar(); |
316 |
|
|
// |
317 |
|
|
// if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on |
318 |
|
|
// |
319 |
|
|
if ( c0->clevel1->good2 == 0 ) { |
320 |
|
|
// |
321 |
|
|
// In selftrigger mode the trkentry variable is set to -1 |
322 |
|
|
// |
323 |
|
|
c0->trkseqno = -3; |
324 |
|
|
// |
325 |
|
|
// Copy values in the class ca from the structure clevel2 |
326 |
|
|
// |
327 |
|
|
c0->FillTrkVar(ca,ntrkentry); |
328 |
|
|
ntrkentry++; |
329 |
|
|
filled = true; |
330 |
|
|
// |
331 |
|
|
} else { |
332 |
|
|
if ( verbose ) printf(" Selftrigger: problems with event at entry %i \n",ye); |
333 |
|
|
}; |
334 |
|
|
// |
335 |
|
|
}; |
336 |
|
|
// |
337 |
|
|
// Call high energy nuclei routine |
338 |
|
|
// |
339 |
|
|
if ( hZn && c0->clevel2->trigty >= 2. ){ |
340 |
|
|
if ( debug ) printf(" Calling selftrigger high energy nuclei routine for entry %i \n",ye); |
341 |
|
|
// |
342 |
|
|
// Disable "track mode" in the fortran routine |
343 |
|
|
// |
344 |
|
|
c0->clevel1->good2 = 0; |
345 |
|
|
// |
346 |
|
|
// Set high energy nuclei flag to one |
347 |
|
|
// |
348 |
|
|
c0->clevel1->hzn = 1; |
349 |
|
|
c0->clevel1->riginput = rigdefault; |
350 |
|
|
// |
351 |
|
|
// We have a selftrigger event to analyze. |
352 |
|
|
// |
353 |
|
|
for (Int_t e = 0; e < 5 ; e++){ |
354 |
|
|
c0->clevel1->al_p[e][0] = 0.; |
355 |
|
|
c0->clevel1->al_p[e][1] = 0.; |
356 |
|
|
}; |
357 |
|
|
c0->clevel1->trkchi2 = 0; |
358 |
|
|
// |
359 |
|
|
c0->GetTrkVar(); |
360 |
|
|
// |
361 |
|
|
// if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on |
362 |
|
|
// |
363 |
|
|
if ( c0->clevel1->good2 == 0 ) { |
364 |
|
|
// |
365 |
|
|
// In selftrigger mode the trkentry variable is set to -1 |
366 |
|
|
// |
367 |
|
|
c0->trkseqno = -2; |
368 |
|
|
// |
369 |
|
|
// Copy values in the class ca from the structure clevel2 |
370 |
|
|
// |
371 |
|
|
c0->FillTrkVar(ca,ntrkentry); |
372 |
|
|
ntrkentry++; |
373 |
|
|
filled = true; |
374 |
|
|
// |
375 |
|
|
} else { |
376 |
|
|
if ( verbose ) printf(" Selftrigger: problems with event at entry %i \n",ye); |
377 |
|
|
}; |
378 |
|
|
// |
379 |
|
|
}; |
380 |
|
|
// |
381 |
|
|
// self trigger event |
382 |
|
|
// |
383 |
|
|
if ( st && c0->clevel2->trigty >= 2. ){ |
384 |
|
|
if ( verbose ) printf(" Selftrigger event at entry %i \n",ye); |
385 |
|
|
// |
386 |
|
|
// Disable "track mode" in the fortran routine |
387 |
|
|
// |
388 |
|
|
c0->clevel1->good2 = 0; |
389 |
|
|
// |
390 |
|
|
// disable high enery nuclei flag; |
391 |
|
|
// |
392 |
|
|
c0->clevel1->hzn = 0; |
393 |
|
|
// |
394 |
|
|
// We have a selftrigger event to analyze. |
395 |
|
|
// |
396 |
|
|
for (Int_t e = 0; e < 5 ; e++){ |
397 |
|
|
c0->clevel1->al_p[e][0] = 0.; |
398 |
|
|
c0->clevel1->al_p[e][1] = 0.; |
399 |
|
|
}; |
400 |
|
|
c0->clevel1->trkchi2 = 0; |
401 |
|
|
// |
402 |
|
|
c0->GetTrkVar(); |
403 |
|
|
// |
404 |
|
|
// if we had no problem (clevel2->good = 0, NOTICE zero, not one in selftrigger mode!), fill and go on |
405 |
|
|
// |
406 |
|
|
if ( c0->clevel1->good2 == 0 ) { |
407 |
|
|
// |
408 |
|
|
// In selftrigger mode the trkentry variable is set to -1 |
409 |
|
|
// |
410 |
|
|
c0->trkseqno = -1; |
411 |
|
|
// |
412 |
|
|
// Copy values in the class ca from the structure clevel2 |
413 |
|
|
// |
414 |
|
|
c0->FillTrkVar(ca,ntrkentry); |
415 |
|
|
ntrkentry++; |
416 |
|
|
filled = true; |
417 |
|
|
// |
418 |
|
|
} else { |
419 |
|
|
if ( verbose ) printf(" Selftrigger: problems with event at entry %i \n",ye); |
420 |
|
|
}; |
421 |
|
|
}; |
422 |
|
|
// |
423 |
|
|
// Clear structures used to communicate with fortran |
424 |
|
|
// |
425 |
|
|
c0->ClearStructs(); |
426 |
|
|
// |
427 |
|
|
// Fill the rootple |
428 |
|
|
// |
429 |
|
|
calo->Fill(); |
430 |
|
|
// |
431 |
|
|
}; |
432 |
|
|
// |
433 |
|
|
}else{ |
434 |
|
|
cout << "Chain entry "<<iev<<" -- ERROR --"<<endl; |
435 |
|
|
}; |
436 |
|
|
}; |
437 |
|
|
// |
438 |
|
|
cout <<endl; |
439 |
|
|
// |
440 |
|
|
event->Clear(); |
441 |
|
|
// |
442 |
|
|
outft->cd(); |
443 |
|
|
event->WriteCloneTrees(); |
444 |
|
|
calo->Write(); |
445 |
|
|
outft->Close(); |
446 |
|
|
// |
447 |
|
|
return; |
448 |
|
|
// |
449 |
|
|
} |