/[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.3 - (hide annotations) (download)
Wed Jul 22 08:18:11 2009 UTC (15 years, 7 months ago) by pizzolot
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +1 -1 lines
adapted to root5.24

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

  ViewVC Help
Powered by ViewVC 1.1.23