/[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.6 - (show annotations) (download)
Mon Mar 12 14:32:46 2007 UTC (17 years, 8 months ago) by pam-rm2
Branch: MAIN
CVS Tags: v2r04, v2r03, HEAD
Changes since 1.5: +18 -17 lines
fixed OBT bug

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

  ViewVC Help
Powered by ViewVC 1.1.23