/[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.1 - (hide annotations) (download)
Thu Jun 15 14:00:30 2006 UTC (18 years, 5 months ago) by pam-rm2
Branch: MAIN
Branch point for: QLflightS4_ND
Initial revision

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     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