/[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.6 - (hide 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 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.6 ULong_t obt;
68 pam-rm2 1.1 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.6 obt=ph->GetOrbitalTime();
142     if(obt<=firstime) firstime=ph->GetOrbitalTime();
143     if(obt>=lastime) lastime=ph->GetOrbitalTime();
144 pam-rm2 1.2 }
145 pam-rm2 1.1
146 pam-rm2 1.6 const ULong_t nint=((lastime-firstime)/(DeltaTevtime*60000));
147     const ULong_t nint1=((lastime-firstime)/(DeltaT));//
148     const ULong_t nint2=lastime-firstime;
149 pam-rm2 1.3 int nbin=(int)nint;
150     int nbin1=(int)nint1;
151     int nbin2=(int)nint2;
152 pam-rm2 1.5 if(nbin2 >= 37000000) nbin2=37000000;
153 pam-rm2 1.6 Double_t obmin=firstime;
154     Double_t obmax=lastime;
155 pam-rm2 1.3
156 pam-rm2 1.1 //************************** HISTOGRAMS ***************************************************//
157     //---------------------------- First histograms -----------------------------------------//
158     oss.str("");
159 pam-rm2 1.2 oss <<"Distribution of Bottom Background Neutrons collected in DeltaT = " << DeltaT << " ms";
160 pam-rm2 1.1 oss1.str("");
161 pam-rm2 1.2 oss1 <<"Distribution of Upper Background Neutrons collected in DeltaT = " << DeltaT << " ms";
162 pam-rm2 1.1 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 pam-rm2 1.2 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 pam-rm2 1.5
173 pam-rm2 1.1 //---------------------------- 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 pam-rm2 1.6 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 pam-rm2 1.5
207 pam-rm2 1.6 Trig->Fill(obt, yTrig);
208     Trigger->Fill(obt, yTrig);
209 pam-rm2 1.1 if(yTrig >=10)
210 pam-rm2 1.6 Triggercut->Fill(obt, yTrig);
211 pam-rm2 1.1 yUpperBackground=0;
212     yBottomBackground=0;
213     yTrig=0;
214     }
215    
216 pam-rm2 1.2
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 pam-rm2 1.1 h1bis->Fill(h1->GetBinContent(i));
227     h3bis->Fill(h3->GetBinContent(i));
228 pam-rm2 1.2 }
229 pam-rm2 1.1
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 pam-rm2 1.5 gPad->SetLogy();
238 pam-rm2 1.1 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 pam-rm2 1.5 gPad->SetLogy();
249 pam-rm2 1.1 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 pam-rm2 1.2 histo1->GetYaxis()->SetTitleOffset(1);
268 pam-rm2 1.1 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 pam-rm2 1.5 histo1->SetMarkerSize(.3);
276 pam-rm2 1.1 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 pam-rm2 1.5 histo2->SetMarkerSize(.3);
284 pam-rm2 1.1 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 pam-rm2 1.5 histo1bis->SetMarkerSize(.3);
297 pam-rm2 1.1 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 pam-rm2 1.5 Trig->SetMarkerSize(.4);
337 pam-rm2 1.1 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 pam-rm2 1.5 Trigger->SetMarkerSize(.4);
353 pam-rm2 1.1 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 pam-rm2 1.5 Triggercut->SetMarkerSize(.4);
370 pam-rm2 1.1 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 pam-rm2 1.2 ULong_t DeltaT =1000;//era 1000
414 pam-rm2 1.4 Double_t DeltaTevtime =1;
415 pam-rm2 1.1
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 pam-rm2 1.2 printf( "\t -DeltaTevtime Time interval to control time evolution of acquired data in minute [default = 5 min.] \n");
428 pam-rm2 1.1 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 pam-rm2 1.2 DeltaTevtime = atof(argv[i]);
473 pam-rm2 1.1 continue;
474     }
475     }
476     }
477    
478     ND_QL(argv[1], outDir, format, DeltaT, DeltaTevtime);
479     }
480    

  ViewVC Help
Powered by ViewVC 1.1.23