/[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.2 - (show annotations) (download)
Fri Jun 23 16:06:09 2006 UTC (18 years, 5 months ago) by pam-rm2
Branch: MAIN
CVS Tags: v1r02, v1r03
Changes since 1.1: +39 -46 lines
corretto il bakground del ND e la rate media di S4

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

  ViewVC Help
Powered by ViewVC 1.1.23