/[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.2 - (hide annotations) (download)
Fri Jun 23 16:06:09 2006 UTC (18 years, 5 months ago) by pam-rm2
Branch: MAIN
CVS Tags: v1r02, v1r03
Changes since 1.1: +39 -46 lines
corretto il bakground del ND e la rate media di S4

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

  ViewVC Help
Powered by ViewVC 1.1.23