source: trunk/MagicSoft/Mars/mtemp/mifae/macros/alpha.C@ 5489

Last change on this file since 5489 was 4117, checked in by rico, 21 years ago
*** empty log message ***
File size: 6.1 KB
Line 
1/* ======================================================================== *\!
2! Author(s): Oscar Blanch & Javier Rico
3\* ======================================================================== */
4
5void alpha(
6 Float_t cutsize=2000,
7 Float_t upcutsize=9999999,
8 Float_t lengthlowcut=0.,
9 Float_t lengthcut=0.26,
10 Float_t widthlowcut=0.,
11 Float_t widthcut=0.12
12)
13{
14 // constants
15 const Int_t nbins = 9;
16 TString psname("alphaMrk.gif");
17 // general settings
18 gROOT->Reset();
19 gStyle->SetCanvasColor(0);
20 gStyle->SetCanvasBorderMode(0);
21 gStyle->SetPadBorderMode(0);
22 gStyle->SetFrameBorderMode(0);
23 gStyle->SetOptTitle(0);
24 gStyle->SetTitleOffset(1.7,"y");
25 gStyle->SetPadLeftMargin(0.15);
26 gStyle->SetOptStat(111110);
27 gStyle->SetStatColor(0);
28 gStyle->SetStatBorderSize(1);
29 gStyle->SetStatW(0.2);
30 gStyle->SetStatH(0.1);
31 gStyle->SetStatX(0.9);
32 gStyle->SetStatY(0.9);
33
34 // define files, chain, trees, canvas...
35 TChain* ton= new TChain("Parameters");
36
37 // ON data
38 ton->Add("/mnt/users/jrico/magic/mars/mars/mtemp/mifae/hillas/crab20040215OnRotateCalA-D.root");
39
40 TChain* toff= new TChain("Parameters");
41
42 // OFF data
43 toff->Add("/mnt/users/jrico/magic/mars/mars/mtemp/mifae/hillas/crab20040215OffRotateCalA-H.root");
44
45 TCanvas *c1 = new TCanvas("c1","c1",800,500);
46 c1->cd(1);
47
48 // define aliases
49 ton->SetAlias("length","MHillas.fLength*0.6/189");
50 ton->SetAlias("width","MHillas.fWidth*0.6/189");
51 ton->SetAlias("dist","MHillasSrc.fDist*0.6/189");
52 ton->SetAlias("conc","MNewImagePar.fConc");
53 ton->SetAlias("size","MHillas.fSize");
54 ton->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
55 ton->SetAlias("alpha","abs(MHillasSrc.fAlpha)");
56 ton->SetAlias("leakage","MNewImagePar.fInnerLeakage1");
57
58 toff->SetAlias("length","MHillas.fLength*0.6/189");
59 toff->SetAlias("width","MHillas.fWidth*0.6/189");
60 toff->SetAlias("dist","MHillasSrc.fDist*0.6/189");
61 toff->SetAlias("conc","MNewImagePar.fConc");
62 toff->SetAlias("size","MHillas.fSize");
63 toff->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
64 toff->SetAlias("alpha","abs(MHillasSrc.fAlpha)");
65 toff->SetAlias("leakage","MNewImagePar.fInnerLeakage1");
66
67 // define cut(s)
68// const char* eventcut="event%4==0 &&";
69 const char* eventcut="";
70 const char varcut[512];
71 // sprintf(varcut,"size>%f && size<%f && length>%f &&length<%f && width>%f && width<%f",cutsize,upcutsize,lengthlowcut,lengthcut,widthlowcut,widthcut);
72 sprintf(varcut,"size>%f && size<%f && dist>0.2 && dist<1.1 && leakage==0 && length>%f &&length<%f && width>%f && width<%f",cutsize,upcutsize,lengthlowcut,lengthcut,widthlowcut,widthcut);
73// sprintf(varcut,"size>%f && dist>0.2 && dist<0.8 && length<0.21 && width<0.125",cutsize,upcutsize,lengthcut,widthcut);
74// sprintf(varcut,"size>%f && size<%f && dist>0.2 && dist<0.8 && length<0.21 && width<0.125",cutsize,upcutsize);
75// sprintf(varcut,"size>%f && size<%f && dist>0.2 && dist<0.8 && length<0.21",cutsize,upcutsize);
76// sprintf(varcut,"size>%f && size<%f && dist>0.2 && dist<0.8",cutsize,upcutsize);
77// sprintf(varcut,"size>%f && size<%f",cutsize,upcutsize)
78// sprintf(varcut,"");
79 Char_t cut[256];
80 Char_t cutoff[256];
81 strcat(cut,eventcut);
82 strcat(cut,varcut);
83 strcat(cutoff,varcut);
84
85 // fill on/off histos
86 TH1F* onhisto = new TH1F("OnHist" ,"Alpha Plot",nbins,0,90);
87 TH1F* offhisto = new TH1F("OffHist","Alpha Plot",nbins,0,90);
88 ton->Draw("alpha>>OnHist",cut);
89 toff->Draw("alpha>>OffHist",cutoff);
90 cout << "Applied cut: " << cut << endl;
91
92 // line/normalization
93 const Int_t inibin = 20./90.*nbins+1;
94 Float_t level=0;
95 Float_t leveloff=0;
96 Float_t levelofferror=0;
97 for(Int_t ibin = inibin; ibin<=nbins;ibin++)
98 {
99 level+=onhisto->GetBinContent(ibin);
100 leveloff+=offhisto->GetBinContent(ibin);
101 }
102 level/=(nbins-inibin+1);
103 leveloff/=(nbins-inibin+1);
104
105 // normalize on/off
106 offhisto->Sumw2(); // needed to compute correct errors after normalization
107 const Float_t norm = level/leveloff;
108 cout << "Normalizing by factor " << norm <<endl;
109 offhisto->Scale(norm);
110
111 // significance:
112 Float_t sig=0,bg=0,esig=0,ebg=0;
113 Float_t significance=0;
114 const Int_t signbins = inibin-1;
115 for(Int_t ibin = 1; ibin<=signbins;ibin++)
116 {
117// Float_t sigma = (onhisto->GetBinContent(ibin)-level)/onhisto->GetBinError(ibin);
118// significance+=sigma*sigma;
119 sig += onhisto->GetBinContent(ibin);
120 esig += onhisto->GetBinError(ibin)*onhisto->GetBinError(ibin);
121 bg += offhisto->GetBinContent(ibin);
122 ebg += offhisto->GetBinError(ibin)*offhisto->GetBinError(ibin);
123 }
124 Float_t error= TMath::Sqrt(esig+ebg);
125 significance = (sig-bg)/error;
126
127 cout << "Excess: " << sig-bg << endl;
128 cout << "N bkg: " << bg << endl;
129 cout << "Significance: " << significance << endl;
130
131 // plot
132 onhisto->SetXTitle("alpha (deg)");
133 onhisto->SetYTitle("Entries");
134 onhisto->SetMarkerColor(4);
135 onhisto->SetMarkerStyle(20);
136 onhisto->SetMarkerSize(.7);
137 onhisto->SetLineColor(4);
138 onhisto->SetMaximum(onhisto->GetBinContent(nbins)*1.5);
139
140 offhisto->SetFillColor(2);
141 offhisto->SetMaximum(offhisto->GetBinContent(nbins)*1.5);
142 offhisto->SetXTitle("alpha (deg)");
143 offhisto->SetYTitle("Entries");
144
145 offhisto->Draw("HIST");
146 offhisto->Draw("ESAME");
147 onhisto->Draw("ESAMES");
148 // onhisto->Draw("E");
149
150 // move stat box to make them all visible
151 gPad->Update();
152 TPaveStats* offpavstat = (TPaveStats*) offhisto->GetListOfFunctions()->FindObject("stats");
153 if(offpavstat)
154 {
155 Float_t shiftx = offpavstat->GetX2NDC()-offpavstat->GetX1NDC();
156 offpavstat->SetX1NDC(offpavstat->GetX1NDC()-shiftx);
157 offpavstat->SetX2NDC(offpavstat->GetX2NDC()-shiftx);
158 }
159
160 // draw line
161 TLine* line = new TLine(0.,level,90.,level);
162 line->SetLineColor(2);
163 line->SetLineStyle(2);
164// line->Draw();
165
166 gPad->Modified();
167 gPad->Update();
168
169 // produce ps
170 // char psname[50];
171 // sprintf(psname,"%s.ps",var);
172 c1->Print(psname);
173}
174
175
Note: See TracBrowser for help on using the repository browser.