/[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.1.1.1 - (show annotations) (download) (vendor branch)
Thu Jun 15 14:00:30 2006 UTC (18 years, 6 months ago) by pam-rm2
Branch: QLflightS4_ND
CVS Tags: v1r01, second
Changes since 1.1: +0 -0 lines
cambiati i grafici con isto, ora lavora su files grandi

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

  ViewVC Help
Powered by ViewVC 1.1.23