/[PAMELA software]/calo/flight/FCaloREPROC/src/CaloREPCore.cpp
ViewVC logotype

Annotation of /calo/flight/FCaloREPROC/src/CaloREPCore.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Fri Oct 26 11:41:32 2007 UTC (17 years, 4 months ago) by mocchiut
Branch: MAIN
Branch point for: FCaloREPROC
Initial revision

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     }

  ViewVC Help
Powered by ViewVC 1.1.23