/[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.5 - (hide annotations) (download)
Mon Sep 25 08:42:56 2006 UTC (18 years, 2 months ago) by pam-rm2
Branch: MAIN
CVS Tags: v2r02, v2r01
Changes since 1.4: +11 -6 lines
ND: introdotta scala logaritmica nel background

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 pam-rm2 1.5 if(nbin2 >= 37000000) nbin2=37000000;
152 pam-rm2 1.3 double obmin=firstime;
153     double obmax=lastime;
154    
155 pam-rm2 1.1 //************************** HISTOGRAMS ***************************************************//
156     //---------------------------- First histograms -----------------------------------------//
157     oss.str("");
158 pam-rm2 1.2 oss <<"Distribution of Bottom Background Neutrons collected in DeltaT = " << DeltaT << " ms";
159 pam-rm2 1.1 oss1.str("");
160 pam-rm2 1.2 oss1 <<"Distribution of Upper Background Neutrons collected in DeltaT = " << DeltaT << " ms";
161 pam-rm2 1.1 TH1F *h1 = new TH1F (filename.Data(), oss.str().c_str(), nbin1,obmin,obmax);
162     TH1F *h3 = new TH1F (filename.Data(), oss1.str().c_str(), nbin1,obmin,obmax);
163     //---------------------------- Second histograms -----------------------------------------//
164 pam-rm2 1.2 title.str("");
165     title<<filename<<": Bottom Background & Upper Background: Time Evolution (DT="<<DeltaTevtime<<" min)";
166     title1.str("");
167     title1<<filename<<". Ratio: UpperBackground / BottomBackground: Time Evolution (DT="<<DeltaTevtime<<" min)";
168     TH1F *histo1= new TH1F(title.str().c_str(),title.str().c_str(),nbin,obmin,obmax);
169     TH1F *histo1bis= new TH1F(title1.str().c_str(),title1.str().c_str(),nbin,obmin,obmax);
170     TH1F *histo2= new TH1F(title.str().c_str(),title.str().c_str(),nbin,obmin,obmax);
171 pam-rm2 1.5
172 pam-rm2 1.1 //---------------------------- Third histograms -----------------------------------------//
173    
174     oss.str("");
175     oss << filename.Data() << ": " << "Trigger events: Time Evolution";
176     oss1.str("");
177     oss1 <<filename.Data() << ": " << "Trigger counting: Time Evolution (DT= "<<DeltaTevtime<<" min)";
178     oss2.str("");
179     oss2 << filename.Data() << ": " << "Trigger counting (events with more than 10 n): Time Evolution (DT= "<<DeltaTevtime<<" min)";
180    
181     TH1F *Trig = new TH1F(oss.str().c_str(), oss.str().c_str(),nbin2,obmin,obmax);
182     TH1F *Trigger = new TH1F(oss1.str().c_str(), oss1.str().c_str(),nbin,obmin,obmax);
183     TH1F *Triggercut = new TH1F(oss2.str().c_str(), oss2.str().c_str(),nbin,obmin,obmax);
184    
185     //-------------------------- Fill histograms from root-ple -----------------------------//
186     for (Int_t i = 0; i < nevents; i++){
187     ndBr->GetEntry(i);
188     headBr->GetEntry(i);
189     tmpSize = ne->Records->GetEntries();
190     if (ne->unpackError == 1) continue; //Check if NeutronEvent is corrupted
191     ph = eh->GetPscuHeader();
192     obt= ph->GetOrbitalTime();
193     for (Int_t j = 0; j < tmpSize; j++){
194     nr = (pamela::neutron::NeutronRecord*)ne->Records->At(j);
195     yTrig = yTrig + (int)nr->trigPhysics;
196     yUpperBackground = yUpperBackground + (int)nr->upperBack;
197     yBottomBackground = yBottomBackground + (int)nr->bottomBack;
198     }
199     histo1->Fill(ph->GetOrbitalTime(), yUpperBackground);
200     histo1bis->Fill(ph->GetOrbitalTime(), yUpperBackground);
201     histo2->Fill(ph->GetOrbitalTime(), yBottomBackground);
202     h1->Fill(ph->GetOrbitalTime(), yBottomBackground);
203     h3->Fill(ph->GetOrbitalTime(), yUpperBackground);
204 pam-rm2 1.5
205 pam-rm2 1.1 Trig->Fill(ph->GetOrbitalTime(), yTrig);
206     Trigger->Fill(ph->GetOrbitalTime(), yTrig);
207     if(yTrig >=10)
208     Triggercut->Fill(ph->GetOrbitalTime(), yTrig);
209     yUpperBackground=0;
210     yBottomBackground=0;
211     yTrig=0;
212     }
213    
214 pam-rm2 1.2
215     Int_t newBin1=(int)((h1->GetMaximum())-(h1->GetMinimum()));
216     Int_t newBin3=(int)((h3->GetMaximum())-(h3->GetMinimum()));
217     oss.str("");
218     oss <<"Distribution of Bottom Background Neutrons collected in DeltaT = " << DeltaT << " ms";
219     oss1.str("");
220     oss1 <<"Distribution of Upper Background Neutrons collected in DeltaT = " << DeltaT << " ms";
221     TH1F *h1bis = new TH1F (filename.Data(), oss.str().c_str(), newBin1, h1->GetMinimum(), h1->GetMaximum());
222     TH1F *h3bis = new TH1F (filename.Data(), oss1.str().c_str(), newBin3, h3->GetMinimum(), h3->GetMaximum());
223     for (Int_t i = 0; i < nbin1; i++){
224 pam-rm2 1.1 h1bis->Fill(h1->GetBinContent(i));
225     h3bis->Fill(h3->GetBinContent(i));
226 pam-rm2 1.2 }
227 pam-rm2 1.1
228     //********************************* CANVASES ************************************************//
229     //--------------------------------- First canvas --------------------------------------------//
230     TCanvas *histoCanv = new TCanvas("ND_QL_1", base.Data(), 1280, 1024);
231     histoCanv->SetFillColor(10);
232     histoCanv->Divide(1,2);
233    
234     histoCanv->cd(1);
235 pam-rm2 1.5 gPad->SetLogy();
236 pam-rm2 1.1 h1bis->SetLineColor(kRed);
237     h1bis->SetFillStyle(3004);
238     h1bis->SetFillColor(kRed);
239     h1bis->GetXaxis()->SetTitle("Number of neutrons");
240     h1bis->GetXaxis()->CenterTitle();
241     h1bis->GetYaxis()->SetTitle("Number of events");
242     h1bis->GetYaxis()->CenterTitle();
243     h1bis->Draw();
244    
245     histoCanv->cd(2);
246 pam-rm2 1.5 gPad->SetLogy();
247 pam-rm2 1.1 h3bis->SetLineColor(kBlue);
248     h3bis->SetFillStyle(3004);
249     h3bis->SetFillColor(kBlue);
250     h3bis->GetXaxis()->SetTitle("Number of neutrons");
251     h3bis->GetXaxis()->CenterTitle();
252     h3bis->GetYaxis()->SetTitle("Number of events");
253     h3bis->GetYaxis()->CenterTitle();
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 pam-rm2 1.2 histo1->GetYaxis()->SetTitleOffset(1);
266 pam-rm2 1.1 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 pam-rm2 1.5 histo1->SetMarkerSize(.3);
274 pam-rm2 1.1 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 pam-rm2 1.5 histo2->SetMarkerSize(.3);
282 pam-rm2 1.1 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 pam-rm2 1.5 histo1bis->SetMarkerSize(.3);
295 pam-rm2 1.1 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 pam-rm2 1.5 Trig->SetMarkerSize(.4);
335 pam-rm2 1.1 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 pam-rm2 1.5 Trigger->SetMarkerSize(.4);
351 pam-rm2 1.1 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 pam-rm2 1.5 Triggercut->SetMarkerSize(.4);
368 pam-rm2 1.1 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 pam-rm2 1.2 ULong_t DeltaT =1000;//era 1000
412 pam-rm2 1.4 Double_t DeltaTevtime =1;
413 pam-rm2 1.1
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 pam-rm2 1.2 printf( "\t -DeltaTevtime Time interval to control time evolution of acquired data in minute [default = 5 min.] \n");
426 pam-rm2 1.1 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 pam-rm2 1.2 DeltaTevtime = atof(argv[i]);
471 pam-rm2 1.1 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