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

Annotation of /quicklook/QLflightS4_ND/ND_QL.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (hide 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 pam-rm2 1.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 pam-rm2 1.2 void ND_QL(TString base, TString outDir, TString format, ULong_t DeltaT, Double_t DeltaTevtime){ // DeltaT in ms, DeltaTevtime in minute
58 pam-rm2 1.1
59     //---------------- Variables initialization --------------------------//
60     Int_t tmpSize=0;
61     Int_t yTrig=0;
62     Int_t yUpperBackground=0;
63     Int_t yBottomBackground=0;
64 pam-rm2 1.4 ULong_t lastime=0, firstime=700000000;
65 pam-rm2 1.1 Double_t scale= 1./DeltaTevtime;
66 pam-rm2 1.2 stringstream title, title1;
67 pam-rm2 1.1 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 pam-rm2 1.3 headBr->GetEntry(0);
135     ph = eh->GetPscuHeader();
136     firstime = ph->GetOrbitalTime();
137    
138     for (Int_t i = 0; i < nevents; i++){
139 pam-rm2 1.2 headBr->GetEntry(i);
140     ph = eh->GetPscuHeader();
141 pam-rm2 1.3 if((ph->GetOrbitalTime()) <= firstime) firstime=ph->GetOrbitalTime();
142     if((ph->GetOrbitalTime()) >= lastime) lastime=ph->GetOrbitalTime();
143 pam-rm2 1.2 }
144 pam-rm2 1.1
145 pam-rm2 1.3 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 pam-rm2 1.1 //************************** HISTOGRAMS ***************************************************//
155     //---------------------------- First histograms -----------------------------------------//
156     oss.str("");
157 pam-rm2 1.2 oss <<"Distribution of Bottom Background Neutrons collected in DeltaT = " << DeltaT << " ms";
158 pam-rm2 1.1 oss1.str("");
159 pam-rm2 1.2 oss1 <<"Distribution of Upper Background Neutrons collected in DeltaT = " << DeltaT << " ms";
160 pam-rm2 1.1 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 pam-rm2 1.2 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 pam-rm2 1.1 //---------------------------- 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 pam-rm2 1.2
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 pam-rm2 1.1 h1bis->Fill(h1->GetBinContent(i));
222     h3bis->Fill(h3->GetBinContent(i));
223 pam-rm2 1.2 }
224 pam-rm2 1.1
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 pam-rm2 1.2 histo1->GetYaxis()->SetTitleOffset(1);
261 pam-rm2 1.1 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 pam-rm2 1.2 ULong_t DeltaT =1000;//era 1000
407 pam-rm2 1.4 Double_t DeltaTevtime =1;
408 pam-rm2 1.1
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 pam-rm2 1.2 printf( "\t -DeltaTevtime Time interval to control time evolution of acquired data in minute [default = 5 min.] \n");
421 pam-rm2 1.1 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 pam-rm2 1.2 DeltaTevtime = atof(argv[i]);
466 pam-rm2 1.1 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