/[PAMELA software]/quicklook/QLflightS4_ND/ND_QL.cpp
ViewVC logotype

Contents of /quicklook/QLflightS4_ND/ND_QL.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (show annotations) (download)
Wed Aug 9 09:56:39 2006 UTC (18 years, 3 months ago) by pam-rm2
Branch: MAIN
CVS Tags: v1r06, v1r05, v2r00
Changes since 1.3: +2 -2 lines
ND: cambiata la reate da 5 a 1 minuto, verificato che non creshi

1 /*****
2 *
3 * Neutron Detector Quick Look
4 * author Marcelli-Malvezzi
5 * Version 1.00 - March 2006
6 *
7 * Not only new release of NeutronDetectorScan.c; all is changed and a new name is requested!
8 *
9 *
10 * Description: The aim of Neutron Detector QL software is to keep under control the time bahaviour of this detector.
11 * It creates three output canvases: the first and the second are relative to Background channels,
12 * while the third one is relative to trigger events ("events" channel) instead.
13 * Each output relates to one PAMELA marshroot. See documentation for a more detailed description of the output.
14 *
15 *
16 * Parameters:
17 * TSTring base - the path to the root directory for the specific Pamela unpack session
18 * There is no default value, without this input the program will not run
19 * TString outDir - the path where to save the output image (Default = ./)
20 * TString format - the format which will be used fo rsave the produced images (Default = "jpg")
21 * Float_t DeltaT - the time interval in ms for ND records colection used for histograms ND_QL_1 plotting (Default = 1000 ms)
22 * Float_t DeltaTevtime - the time interval in minute for calculation of average ND count rate per minute,
23 * see ND_QL_2 and ND_QL_3 plots (Default = 1 minute)
24 * Version 1.1 - June 2006
25 *
26 * Fixed bugs: for a large namber of events is not possible to initialize vectors, so all graphs have been converted in histograms
27 *
28 */
29
30 #include <iostream>
31 #include <string>
32 #include <sstream>
33 #include <math.h>
34 #include "TStyle.h"
35 #include "TLatex.h"
36 #include "TFile.h"
37 #include "TList.h"
38 #include "TTree.h"
39 #include "TObjString.h"
40 #include "TCanvas.h"
41 #include "TGraph.h"
42 #include "TLegend.h"
43 #include "TH1F.h"
44 #include "TF1.h"
45 #include "TGaxis.h"
46 #include "TString.h"
47 #include "TPaveText.h"
48 #include "EventHeader.h"
49 #include "PscuHeader.h"
50 #include "TMultiGraph.h"
51 #include "physics/neutronDetector/NeutronEvent.h"
52 #include "physics/neutronDetector/NeutronRecord.h"
53
54
55 using namespace std;
56
57 void ND_QL(TString base, TString outDir, TString format, ULong_t DeltaT, Double_t DeltaTevtime){ // DeltaT in ms, DeltaTevtime in minute
58
59 //---------------- Variables initialization --------------------------//
60 Int_t tmpSize=0;
61 Int_t yTrig=0;
62 Int_t yUpperBackground=0;
63 Int_t yBottomBackground=0;
64 ULong_t lastime=0, firstime=700000000;
65 Double_t scale= 1./DeltaTevtime;
66 stringstream title, title1;
67 double obt;
68 stringstream oss, noentries;
69 stringstream oss1, oss2, oss3;
70
71 //------- load root file --------------
72 TFile *file = new TFile(base.Data()) ;
73
74 TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();
75 filename = ((TObjString*)filename.Tokenize('.')->First())->GetString();
76
77 if ( outDir == "" ) outDir = ".";
78 if (!file){
79 printf("No such file in the directory has been found");
80 return;
81 }
82
83 //------------------- Takes the tree of ND ---------------------------//
84 TTree *phystr = (TTree*)file->Get("Physics");
85 TBranch *ndBr = phystr->GetBranch("Neutron");
86 TBranch *headBr = phystr->GetBranch("Header");
87
88 //PAMELA classes
89 pamela::neutron::NeutronEvent *ne=0;
90 pamela::neutron::NeutronRecord *nr=0;
91 pamela::EventHeader *eh=0;
92 pamela::PscuHeader *ph=0;
93
94 phystr->SetBranchAddress("Neutron", &ne);
95 phystr->SetBranchAddress("Header", &eh);
96
97 Long64_t nevents = phystr->GetEntries();
98 const Int_t size = nevents;
99
100 //--------------- output if there are 0 or <0 entries ---------------------------------------//
101 if (nevents<=0){
102 cout<<"WARNING: No Entries, RETURN"<<"\n";
103
104 TCanvas *canvas4 = new TCanvas("No entries", "No entries ", 400, 200);
105 canvas4->SetFillColor(10);
106 canvas4->cd();
107
108 TLatex *l = new TLatex();
109 l->SetTextAlign(12);
110 l->SetTextSize(0.15);
111 l->SetTextColor(2);
112 noentries.str("");
113 noentries<< "ND_QL:";
114 l->DrawLatex(0.05, 0.7, noentries.str().c_str());
115 noentries.str("");
116 noentries<< "No entries for this files";
117 l->DrawLatex(0.05, 0.5, noentries.str().c_str());
118
119 if (outDir == "./") {
120 oss.str("");
121 oss << filename.Data() << "ND_QL." << format.Data();
122 } else {
123 oss.str("");
124 oss << outDir.Data() << filename.Data() << "_ND_QL." << format.Data();
125 }
126
127 canvas4->Update();
128 canvas4->SaveAs(oss.str().c_str());
129
130 return;
131 }
132
133 //------------ first and last events -> nbin ---------------------------------//
134 headBr->GetEntry(0);
135 ph = eh->GetPscuHeader();
136 firstime = ph->GetOrbitalTime();
137
138 for (Int_t i = 0; i < nevents; i++){
139 headBr->GetEntry(i);
140 ph = eh->GetPscuHeader();
141 if((ph->GetOrbitalTime()) <= firstime) firstime=ph->GetOrbitalTime();
142 if((ph->GetOrbitalTime()) >= lastime) lastime=ph->GetOrbitalTime();
143 }
144
145 const Double_t nint=((lastime-firstime)/(DeltaTevtime*60000));
146 const Double_t nint1=((lastime-firstime)/(DeltaT));//
147 const Double_t nint2=lastime-firstime;
148 int nbin=(int)nint;
149 int nbin1=(int)nint1;
150 int nbin2=(int)nint2;
151 double obmin=firstime;
152 double obmax=lastime;
153
154 //************************** HISTOGRAMS ***************************************************//
155 //---------------------------- First histograms -----------------------------------------//
156 oss.str("");
157 oss <<"Distribution of Bottom Background Neutrons collected in DeltaT = " << DeltaT << " ms";
158 oss1.str("");
159 oss1 <<"Distribution of Upper Background Neutrons collected in DeltaT = " << DeltaT << " ms";
160 TH1F *h1 = new TH1F (filename.Data(), oss.str().c_str(), nbin1,obmin,obmax);
161 TH1F *h3 = new TH1F (filename.Data(), oss1.str().c_str(), nbin1,obmin,obmax);
162 //---------------------------- Second histograms -----------------------------------------//
163 title.str("");
164 title<<filename<<": Bottom Background & Upper Background: Time Evolution (DT="<<DeltaTevtime<<" min)";
165 title1.str("");
166 title1<<filename<<". Ratio: UpperBackground / BottomBackground: Time Evolution (DT="<<DeltaTevtime<<" min)";
167 TH1F *histo1= new TH1F(title.str().c_str(),title.str().c_str(),nbin,obmin,obmax);
168 TH1F *histo1bis= new TH1F(title1.str().c_str(),title1.str().c_str(),nbin,obmin,obmax);
169 TH1F *histo2= new TH1F(title.str().c_str(),title.str().c_str(),nbin,obmin,obmax);
170 //---------------------------- Third histograms -----------------------------------------//
171
172 oss.str("");
173 oss << filename.Data() << ": " << "Trigger events: Time Evolution";
174 oss1.str("");
175 oss1 <<filename.Data() << ": " << "Trigger counting: Time Evolution (DT= "<<DeltaTevtime<<" min)";
176 oss2.str("");
177 oss2 << filename.Data() << ": " << "Trigger counting (events with more than 10 n): Time Evolution (DT= "<<DeltaTevtime<<" min)";
178
179 TH1F *Trig = new TH1F(oss.str().c_str(), oss.str().c_str(),nbin2,obmin,obmax);
180 TH1F *Trigger = new TH1F(oss1.str().c_str(), oss1.str().c_str(),nbin,obmin,obmax);
181 TH1F *Triggercut = new TH1F(oss2.str().c_str(), oss2.str().c_str(),nbin,obmin,obmax);
182
183 //-------------------------- Fill histograms from root-ple -----------------------------//
184 for (Int_t i = 0; i < nevents; i++){
185 ndBr->GetEntry(i);
186 headBr->GetEntry(i);
187 tmpSize = ne->Records->GetEntries();
188 if (ne->unpackError == 1) continue; //Check if NeutronEvent is corrupted
189 ph = eh->GetPscuHeader();
190 obt= ph->GetOrbitalTime();
191 for (Int_t j = 0; j < tmpSize; j++){
192 nr = (pamela::neutron::NeutronRecord*)ne->Records->At(j);
193 yTrig = yTrig + (int)nr->trigPhysics;
194 yUpperBackground = yUpperBackground + (int)nr->upperBack;
195 yBottomBackground = yBottomBackground + (int)nr->bottomBack;
196 }
197 histo1->Fill(ph->GetOrbitalTime(), yUpperBackground);
198 histo1bis->Fill(ph->GetOrbitalTime(), yUpperBackground);
199 histo2->Fill(ph->GetOrbitalTime(), yBottomBackground);
200 h1->Fill(ph->GetOrbitalTime(), yBottomBackground);
201 h3->Fill(ph->GetOrbitalTime(), yUpperBackground);
202 Trig->Fill(ph->GetOrbitalTime(), yTrig);
203 Trigger->Fill(ph->GetOrbitalTime(), yTrig);
204 if(yTrig >=10)
205 Triggercut->Fill(ph->GetOrbitalTime(), yTrig);
206 yUpperBackground=0;
207 yBottomBackground=0;
208 yTrig=0;
209 }
210
211
212 Int_t newBin1=(int)((h1->GetMaximum())-(h1->GetMinimum()));
213 Int_t newBin3=(int)((h3->GetMaximum())-(h3->GetMinimum()));
214 oss.str("");
215 oss <<"Distribution of Bottom Background Neutrons collected in DeltaT = " << DeltaT << " ms";
216 oss1.str("");
217 oss1 <<"Distribution of Upper Background Neutrons collected in DeltaT = " << DeltaT << " ms";
218 TH1F *h1bis = new TH1F (filename.Data(), oss.str().c_str(), newBin1, h1->GetMinimum(), h1->GetMaximum());
219 TH1F *h3bis = new TH1F (filename.Data(), oss1.str().c_str(), newBin3, h3->GetMinimum(), h3->GetMaximum());
220 for (Int_t i = 0; i < nbin1; i++){
221 h1bis->Fill(h1->GetBinContent(i));
222 h3bis->Fill(h3->GetBinContent(i));
223 }
224
225 //********************************* CANVASES ************************************************//
226 //--------------------------------- First canvas --------------------------------------------//
227 TCanvas *histoCanv = new TCanvas("ND_QL_1", base.Data(), 1280, 1024);
228 histoCanv->SetFillColor(10);
229 histoCanv->Divide(1,2);
230
231 histoCanv->cd(1);
232 h1bis->SetLineColor(kRed);
233 h1bis->SetFillStyle(3004);
234 h1bis->SetFillColor(kRed);
235 h1bis->GetXaxis()->SetTitle("Number of neutrons");
236 h1bis->GetXaxis()->CenterTitle();
237 h1bis->GetYaxis()->SetTitle("Number of events");
238 h1bis->GetYaxis()->CenterTitle();
239 h1bis->Draw();
240
241 histoCanv->cd(2);
242 h3bis->SetLineColor(kBlue);
243 h3bis->SetFillStyle(3004);
244 h3bis->SetFillColor(kBlue);
245 h3bis->GetXaxis()->SetTitle("Number of neutrons");
246 h3bis->GetXaxis()->CenterTitle();
247 h3bis->GetYaxis()->SetTitle("Number of events");
248 h3bis->GetYaxis()->CenterTitle();
249 h3bis->Draw();
250
251 //---------------------------------- Second Canvas -----------------------------------------//
252 TCanvas *finalCanv = new TCanvas("ND_QL_2", base.Data(), 1280, 1024);
253 finalCanv->SetFillColor(10);
254 finalCanv->Divide(1,2);
255
256 finalCanv->cd(1);
257 histo1->SetStats(kFALSE);
258 histo1->GetYaxis()->SetTitle("Number of recorded Neutrons / DT (min^-1) ");
259 histo1->GetYaxis()->SetTitleSize(0.04);
260 histo1->GetYaxis()->SetTitleOffset(1);
261 histo1->GetYaxis()->CenterTitle();
262 histo1->GetYaxis()->SetLabelSize(0.03);
263 histo1->GetXaxis()->CenterTitle();
264 histo1->GetXaxis()->SetTitleSize(0.04);
265 histo1->GetXaxis()->SetTitleOffset(1);
266 histo1->GetXaxis()->SetLabelSize(0.03);
267 histo1->GetXaxis()->SetTitle("OBT (ms)");
268 histo1->SetMarkerSize(.5);
269 histo1->SetMarkerStyle(21);
270 histo1->SetMarkerColor(3);
271 histo1->Scale(scale);
272 histo1->SetMinimum(0.8);
273 histo1->Draw("9p0");
274
275 histo2->SetStats(kFALSE);
276 histo2->SetMarkerSize(.5);
277 histo2->SetMarkerStyle(21);
278 histo2->SetMarkerColor(2);
279 histo2->Scale(scale);
280 histo2->SetMinimum(0.8);
281 histo2->Draw("9p0same");
282
283 TLegend *leg = new TLegend(.91,.65, .99, .81);
284 leg->AddEntry(histo2,"Bottom");
285 leg->AddEntry(histo1,"Upper");
286 leg->Draw();
287
288 finalCanv->cd(2);
289 histo1bis->SetMarkerSize(.5);
290 histo1bis->SetStats(kFALSE);
291 histo1bis->SetMarkerStyle(21);
292 histo1bis->SetMarkerColor(4);
293 histo1bis->GetXaxis()->SetTitle("OBT (ms)");
294 histo1bis->GetXaxis()->SetTitleSize(0.04);
295 histo1bis->GetXaxis()->CenterTitle();
296 histo1bis->GetXaxis()->SetLabelSize(0.03);
297 histo1bis->GetYaxis()->SetTitle("Upper / Bottom Background");
298 histo1bis->GetYaxis()->CenterTitle();
299 histo1bis->GetYaxis()->SetTitleSize(0.04);
300 histo1bis->GetYaxis()->SetTitleOffset(0.5);
301 histo1bis->GetYaxis()->SetLabelSize(0.03);
302 histo1bis->Scale(scale);
303 histo1bis->Divide(histo2);
304 histo1bis->SetMinimum(0.);
305 if(histo1bis->GetMaximum()<2)
306 histo1bis->SetMaximum(2.0);
307 histo1bis->Draw("9p0");
308 TF1 *func1 = new TF1("func1", "1.4");
309 func1->SetRange(firstime, lastime);
310 func1->SetLineColor(3);
311 func1->SetLineStyle(1);
312 func1->SetLineWidth(3);
313 func1->Draw("same");
314 TF1 *func2 = new TF1("func2", "0.6");
315 func2->SetRange(firstime, lastime);
316 func2->SetLineColor(3);
317 func2->SetLineStyle(1);
318 func2->SetLineWidth(3);
319 func2->Draw("same");
320
321 //---------------------------------- Third Canvas -----------------------------------------//
322 TCanvas *triggerCanv = new TCanvas("ND_QL_3", base.Data(), 1280, 1024);
323 triggerCanv->SetFillColor(10);
324 triggerCanv->Divide(1,3);
325
326 triggerCanv->cd(1);
327 Trig->SetStats(kFALSE);
328 Trig->SetMarkerStyle(21);
329 Trig->SetMarkerSize(.7);
330 Trig->SetMarkerColor(2);
331 Trig->GetXaxis()->SetTitle("OBT (ms)");
332 Trig->GetXaxis()->CenterTitle();
333 Trig->GetXaxis()->SetTitleSize(0.05);
334 Trig->GetXaxis()->SetLabelSize(0.05);
335 Trig->GetYaxis()->SetTitle("Number of Neutrons");
336 Trig->GetYaxis()->CenterTitle();
337 Trig->GetYaxis()->SetTitleSize(0.06);
338 Trig->GetYaxis()->SetTitleOffset(0.5);
339 Trig->GetYaxis()->SetLabelSize(0.05);
340 Trig->Draw("9p0");
341
342 triggerCanv->cd(2);
343 Trigger->SetStats(kFALSE);
344 Trigger->SetMarkerStyle(21);
345 Trigger->SetMarkerSize(.7);
346 Trigger->SetMarkerColor(4);
347 Trigger->GetYaxis()->SetTitle("Number of Neutrons");
348 Trigger->GetYaxis()->SetTitleSize(0.06);
349 Trigger->GetYaxis()->SetTitleOffset(0.5);
350 Trigger->GetYaxis()->CenterTitle();
351 Trigger->GetYaxis()->SetLabelSize(0.05);
352 Trigger->GetXaxis()->CenterTitle();
353 Trigger->GetXaxis()->SetTitleSize(0.05);
354 Trigger->GetXaxis()->SetLabelSize(0.05);
355 Trigger->GetXaxis()->SetTitle("OBT (ms)");
356 Trigger->SetMinimum(0.);
357 Trigger->Draw("9p0");
358
359 triggerCanv->cd(3);
360 Triggercut->SetStats(kFALSE);
361 Triggercut->SetMarkerStyle(21);
362 Triggercut->SetMarkerSize(.7);
363 Triggercut->SetMarkerColor(4);
364 Triggercut->GetYaxis()->SetTitle("Number of Neutrons");
365 Triggercut->GetYaxis()->SetTitleSize(0.06);
366 Triggercut->GetYaxis()->SetTitleOffset(.5);
367 Triggercut->GetYaxis()->CenterTitle();
368 Triggercut->GetYaxis()->SetLabelSize(0.05);
369 Triggercut->GetXaxis()->CenterTitle();
370 Triggercut->GetXaxis()->SetTitleSize(0.05);
371 Triggercut->GetXaxis()->SetTitleOffset(1);
372 Triggercut->GetXaxis()->SetLabelSize(0.05);
373 Triggercut->GetXaxis()->SetTitle("OBT (ms)");
374 Triggercut->Draw("9p");
375
376 //************************** to save canvas ******************************************//
377
378 if (outDir == "") {
379 oss1.str("");
380 oss2.str("");
381 oss3.str("");
382 oss1 << filename.Data() << "." << format.Data();
383 oss2 << filename.Data() << "." << format.Data();
384 oss3 << filename.Data() << "." << format.Data();
385 } else {
386 oss1.str("");
387 oss2.str("");
388 oss3.str("");
389 oss1 << outDir.Data() << filename.Data() << "_ND_QL_1." << format.Data();
390 oss2 << outDir.Data() << filename.Data() << "_ND_QL_2." << format.Data();
391 oss3 << outDir.Data() << filename.Data() << "_ND_QL_3." << format.Data();
392 }
393
394 finalCanv->SaveAs(oss2.str().c_str());
395 histoCanv->SaveAs(oss1.str().c_str());
396 triggerCanv->SaveAs(oss3.str().c_str());
397
398 file->Close();
399
400 }
401
402 int main(int argc, char* argv[]){
403 TString path;
404 TString outDir ="./";
405 TString format ="jpg";
406 ULong_t DeltaT =1000;//era 1000
407 Double_t DeltaTevtime =1;
408
409 if (argc < 2){
410 printf("You have to insert at least the file to analyze \n");
411 printf("Try '--help' for more information. \n");
412 exit(1);
413 }
414 if (!strcmp(argv[1], "--help")){
415 printf( "Usage: ND_QL FILE [OPTION] \n");
416 printf( "\t --help Print this help and exit \n");
417 printf( "\t -outDir[path] Path where to put the output [default ./] \n");
418 printf( "\t -format[jpg|gif|ps] Format for output files [default 'jpg'] \n");
419 printf( "\t -DeltaT Time interval to bin histograms in ms [default = 1000ms] \n");
420 printf( "\t -DeltaTevtime Time interval to control time evolution of acquired data in minute [default = 5 min.] \n");
421 exit(1);
422 }
423 path=argv[1];
424 for (int i = 2; i < argc; i++){
425 if (!strcmp(argv[i], "-outDir")){
426 if (++i >= argc){
427 printf( "-outDir needs arguments. \n");
428 printf( "Try '--help' for more information. \n");
429 exit(1);
430 }
431 else{
432 outDir = argv[i];
433 continue;
434 }
435 }
436 if (!strcmp(argv[i], "-format")){
437 if (++i >= argc){
438 printf( "-format needs arguments. \n");
439 printf( "Try '--help' for more information. \n");
440 exit(1);
441 }
442 else{
443 format = argv[i];
444 continue;
445 }
446 }
447 if (!strcmp(argv[i], "-DeltaT")){
448 if (++i >= argc){
449 printf( "-DeltaT needs arguments. \n");
450 printf( "Try '--help' for more information. \n");
451 exit(1);
452 }
453 else{
454 DeltaT = atol(argv[i]);
455 continue;
456 }
457 }
458 if (!strcmp(argv[i], "-DeltaTevtime")){
459 if (++i >= argc){
460 printf( "-DeltaTevtime needs arguments. \n");
461 printf( "Try '--help' for more information. \n");
462 exit(1);
463 }
464 else{
465 DeltaTevtime = atof(argv[i]);
466 continue;
467 }
468 }
469 }
470
471 ND_QL(argv[1], outDir, format, DeltaT, DeltaTevtime);
472 }
473
474

  ViewVC Help
Powered by ViewVC 1.1.23