/[PAMELA software]/calo/flight/CaloProfile/src/CaloProfile.cpp
ViewVC logotype

Annotation of /calo/flight/CaloProfile/src/CaloProfile.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Fri Nov 9 09:11:24 2007 UTC (17 years, 2 months ago) by mocchiut
Branch: MAIN
Branch point for: v1r00
Initial revision

1 mocchiut 1.1 #include <CaloProfile.h>
2    
3     ////////////////////////////////////////////////////////////////////////
4     /**
5     * 1-dimension function describing lateral distribution of the
6     * shower as viewed by calorimeter
7     * (projection of 3d function in one direction)
8     *
9     * xi[0] = x or y coordinate relative to shower axis
10     * parmin[0] = rt
11     * parmin[1] = p
12     * parmin[2] = rc
13     *
14     */
15     ////////////////////////////////////////////////////////////////////////
16     Double_t cfradx(Double_t *xi, Double_t *parmin) {
17    
18     double fradxmin2,p,rt,rc,es,x,pig,norm,c;
19     x=*xi;
20     pig = acos(-1.);
21     rt=parmin[0];
22     p=parmin[1];
23     rc=parmin[2];
24     norm=parmin[3];
25     c=parmin[4];
26     x=x-c;
27     es=1.5;
28     fradxmin2=p*pig*pow(rc,2)/pow((pow(x,2)+pow(rc,2)),es);
29     fradxmin2=fradxmin2+(1-p)*pig*pow(rt,2)/pow((pow(x,2)+pow(rt,2)),es);
30     fradxmin2=norm*fradxmin2/(2*pig);
31     //cout<<"x,fradxmin2 "<< x<<" "<<fradxmin2 <<endl;
32     return fradxmin2;
33     }
34    
35    
36     Double_t cfunc(Double_t *tin,Double_t *par)
37     {
38    
39     Double_t gaa,a,tm,et,value,t,t0;
40     t=*tin;
41     et=par[0];
42     a=par[1];
43     tm=par[2];
44     t0=par[3];
45     gaa=TMath::Gamma(a);
46    
47     value=et*((a-1)/(tm*gaa))*pow(((a-1)*(t-t0)/tm),(a-1))*exp(-(a-1)*(t-t0)/tm);
48     return value;
49     }
50    
51    
52     //--------------------------------------
53     /**
54     * Default constructor
55     */
56     CaloLat::CaloLat(){
57     Clear();
58     };
59    
60     /**
61     * Default constructor
62     */
63     CaloLong::CaloLong(){
64     Clear();
65     };
66    
67     CaloLat::CaloLat(PamLevel2 *l2p){
68     //
69     Clear();
70     //
71     L2 = l2p;
72     //
73     if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
74     //
75     OBT = 0;
76     PKT = 0;
77     atime = 0;
78     //
79     debug = false;
80     //
81     };
82    
83     CaloLong::CaloLong(PamLevel2 *l2p){
84     //
85     Clear();
86     //
87     L2 = l2p;
88     //
89     if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
90     //
91     OBT = 0;
92     PKT = 0;
93     atime = 0;
94     //
95     debug = false;
96     //
97     };
98    
99     void CaloLat::Clear(){
100     //
101     };
102    
103     void CaloLong::Clear(){
104     //
105     };
106    
107     void CaloLat::Print(){
108     //
109     Process();
110     //
111     printf("==================== Calorimeter Lateral Profile =======================\n");
112     printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
113     // printf(" nx [number of X combination]:.. %i\n",nx);
114     // printf(" ny [number of Y combination]:.. %i\n",ny);
115     printf("========================================================================\n");
116     //
117     };
118    
119     void CaloLong::Print(){
120     //
121     Process();
122     //
123     printf("==================== Calorimeter Lateral Profile =======================\n");
124     printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
125     // printf(" nx [number of X combination]:.. %i\n",nx);
126     // printf(" ny [number of Y combination]:.. %i\n",ny);
127     printf("========================================================================\n");
128     //
129     };
130    
131     void CaloLat::Draw(){
132     //
133     Process();
134     Draw(-1,-1);
135     };
136    
137     void CaloLong::Draw(){
138     //
139     Process();
140     Draw(-1);
141     };
142    
143    
144     void CaloLat::Draw(Int_t view,Int_t plane){
145     //
146     Int_t minv = 0;
147     Int_t maxv = 0;
148     Int_t minp = 0;
149     Int_t maxp = 0;
150     //
151     if ( view == -1 ){
152     minv = 0;
153     maxv = 2;
154     } else {
155     minv = view;
156     maxv = view+1;
157     };
158    
159     if ( plane == -1 ){
160     minp = 0;
161     maxp = 22;
162     } else {
163     minp = plane;
164     maxp = plane+1;
165     };
166     //
167     Process();
168     //
169     gStyle->SetLabelSize(0.04);
170     gStyle->SetNdivisions(510,"XY");
171     //
172     for (Int_t v=minv; v<maxv;v++){
173     for (Int_t p=minp; p<maxp;p++){
174     TString hid = Form("clatv%ip%i",v,p);
175     TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
176     if ( tc ){
177     // tc->Clear();
178     } else {
179     tc = new TCanvas(hid,hid);
180     };
181     //
182     TString thid = Form("hlatv%ip%i",v,p);
183     TH1F *th = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
184     if ( th ) th->Delete();
185     // th->Clear();
186     // th->Reset();
187     // } else {
188     th = new TH1F(thid,thid,96,-0.5,95.5);
189     // };
190     tc->cd();
191     //
192     for (Int_t st=0;st<96;st++){
193     th->Fill(st,estrip[v][p][st]);
194     };
195     th->Draw();
196     tc->Modified();
197     tc->Update();
198     };
199     };
200     //
201     gStyle->SetLabelSize(0);
202     gStyle->SetNdivisions(1,"XY");
203     //
204     };
205    
206     void CaloLong::Draw(Int_t view){
207     //
208     Int_t minv = 0;
209     Int_t maxv = 0;
210     //
211     if ( view == -1 ){
212     maxv = -1;
213     } else {
214     minv = view;
215     maxv = view+1;
216     };
217     //
218     Process();
219     //
220     gStyle->SetLabelSize(0.04);
221     gStyle->SetNdivisions(510,"XY");
222     //
223     if ( maxv != -1 ){
224     for (Int_t v=minv; v<maxv;v++){
225     TString hid = Form("clongv%i",v);
226     TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
227     if ( tc ){
228     // tc->Clear();
229     } else {
230     tc = new TCanvas(hid,hid);
231     };
232     //
233     TString thid = Form("hlongv%i",v);
234     TH1F *th = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
235     if ( th ) th->Delete();
236     // th->Clear();
237     // th->Reset();
238     // } else {
239     th = new TH1F(thid,thid,22,-0.5,21.5);
240     // };
241     tc->cd();
242     //
243     for (Int_t st=0;st<22;st++){
244     th->Fill(st,eplane[v][st]);
245     };
246     th->Draw();
247     tc->Modified();
248     tc->Update();
249     };
250     } else {
251     //
252     TString hid = Form("clongvyvx");
253     TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
254     if ( tc ){
255     } else {
256     tc = new TCanvas(hid,hid);
257     };
258     //
259     TString thid = Form("hlongvyvx");
260     TH1F *th = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
261     if ( th ) th->Delete();
262     th = new TH1F(thid,thid,44,-0.5,43.5);
263     tc->cd();
264     Int_t pp=0;
265     for (Int_t st=0;st<22;st++){
266     for (Int_t v=1; v>=0;v--){
267     //
268     th->Fill(pp,eplane[v][st]);
269     //
270     pp++;
271     };
272     };
273     th->Draw();
274     tc->Modified();
275     tc->Update();
276     };
277     //
278     gStyle->SetLabelSize(0);
279     gStyle->SetNdivisions(1,"XY");
280     //
281     };
282    
283     void CaloLat::Delete(){
284     Clear();
285     //delete this;
286     };
287    
288     void CaloLong::Delete(){
289     Clear();
290     //delete this;
291     };
292    
293     void CaloLat::Process(){
294     //
295     if ( !L2 ){
296     printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
297     printf(" ERROR: CaloHough variables not filled \n");
298     return;
299     };
300     //
301     Bool_t newentry = false;
302     //
303     if ( L2->IsORB() ){
304     if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
305     newentry = true;
306     OBT = L2->GetOrbitalInfo()->OBT;
307     PKT = L2->GetOrbitalInfo()->pkt_num;
308     atime = L2->GetOrbitalInfo()->absTime;
309     };
310     } else {
311     newentry = true;
312     };
313     //
314     if ( !newentry ) return;
315     //
316     if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
317     //
318     Clear();
319     //
320     // let's start
321     //
322     memset(estrip,0, 4224*sizeof(Float_t));
323     Float_t mip1 = 0.;
324     Int_t view1 = 0;
325     Int_t plane1 = 0;
326     Int_t strip1 = 0;
327     //
328     for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
329     mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
330     estrip[view1][plane1][strip1] = mip1;
331     };
332     //
333     if ( debug ) this->Print();
334     if ( debug ) printf(" exit \n");
335     //
336     };
337    
338     void CaloLong::Process(){
339     //
340     if ( !L2 ){
341     printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
342     printf(" ERROR: CaloHough variables not filled \n");
343     return;
344     };
345     //
346     Bool_t newentry = false;
347     //
348     if ( L2->IsORB() ){
349     if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
350     newentry = true;
351     OBT = L2->GetOrbitalInfo()->OBT;
352     PKT = L2->GetOrbitalInfo()->pkt_num;
353     atime = L2->GetOrbitalInfo()->absTime;
354     };
355     } else {
356     newentry = true;
357     };
358     //
359     if ( !newentry ) return;
360     //
361     if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
362     //
363     Clear();
364     //
365     // let's start
366     //
367     memset(eplane,0, 2*22*sizeof(Float_t));
368     //
369     for (Int_t v=0; v<2; v++){
370     for (Int_t i=0; i<22; i++){
371     eplane[v][i] = L2->GetCaloLevel1()->qtotpl(v,i);
372     };
373     };
374     //
375     if ( debug ) this->Print();
376     if ( debug ) printf(" exit \n");
377     //
378     };
379    

  ViewVC Help
Powered by ViewVC 1.1.23