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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.8 - (show annotations) (download)
Tue Aug 11 14:23:09 2009 UTC (15 years, 6 months ago) by mocchiut
Branch: MAIN
Changes since 1.7: +114 -4 lines
New small features added

1 #include <CaloProfile.h>
2 //
3 ClassImp(CaloLat);
4 ClassImp(CaloLong);
5 ClassImp(Calo2D);
6
7 ////////////////////////////////////////////////////////////////////////
8 /**
9 * 1-dimension function describing lateral distribution of the
10 * shower as viewed by calorimeter
11 * (projection of 3d function in one direction)
12 *
13 * xi[0] = x or y coordinate relative to shower axis
14 * parmin[0] = rt
15 * parmin[1] = p
16 * parmin[2] = rc
17 *
18 */
19 ////////////////////////////////////////////////////////////////////////
20 Double_t cfradx(Double_t *xi, Double_t *parmin) {
21
22 double fradxmin2,p,rt,rc,es,x,pig,norm,c;
23 x=*xi;
24 pig = acos(-1.);
25 rt=parmin[0];
26 p=parmin[1];
27 rc=parmin[2];
28 norm=parmin[3];
29 c=parmin[4];
30 x=x-c;
31 es=1.5;
32 fradxmin2=p*pig*pow(rc,2)/pow((pow(x,2)+pow(rc,2)),es);
33 fradxmin2=fradxmin2+(1-p)*pig*pow(rt,2)/pow((pow(x,2)+pow(rt,2)),es);
34 fradxmin2=norm*fradxmin2/(2*pig);
35 //cout<<"x,fradxmin2 "<< x<<" "<<fradxmin2 <<endl;
36 return fradxmin2;
37 }
38
39
40 Double_t cfunc(Double_t *tin,Double_t *par)
41 {
42
43 Double_t gaa,a,tm,et,value,t,t0;
44 t=*tin;
45 et=par[0];
46 a=par[1];
47 tm=par[2];
48 t0=par[3];
49 gaa=TMath::Gamma(a);
50
51 value=et*((a-1)/(tm*gaa))*pow(((a-1)*(t-t0)/tm),(a-1))*exp(-(a-1)*(t-t0)/tm);
52 return value;
53 }
54
55
56 //--------------------------------------
57 /**
58 * Default constructor
59 */
60 CaloLat::CaloLat(){
61 Clear();
62 };
63
64 /**
65 * Default constructor
66 */
67 Calo2D::Calo2D(){
68 Clear();
69 };
70
71 CaloLat::CaloLat(PamLevel2 *l2p){
72 //
73 Clear();
74 //
75 L2 = l2p;
76 //
77 if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
78 //
79 OBT = 0;
80 PKT = 0;
81 atime = 0;
82 //
83 debug = false;
84 //
85 };
86
87 Calo2D::Calo2D(PamLevel2 *l2p){
88 //
89 Clear();
90 //
91 L2 = l2p;
92 //
93 if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
94 //
95 OBT = 0;
96 PKT = 0;
97 atime = 0;
98 //
99 debug = false;
100 //
101 };
102
103 void CaloLat::Clear(){
104 //
105 };
106
107 void Calo2D::Clear(){
108 //
109 };
110
111 void CaloLat::Print(){
112 //
113 Process();
114 //
115 printf("==================== Calorimeter Lateral Profile =======================\n");
116 printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
117 // printf(" nx [number of X combination]:.. %i\n",nx);
118 // printf(" ny [number of Y combination]:.. %i\n",ny);
119 printf("========================================================================\n");
120 //
121 };
122
123 void Calo2D::Print(){
124 //
125 Process();
126 //
127 printf("==================== Calorimeter 2D Profile =======================\n");
128 printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
129 // printf(" nx [number of X combination]:.. %i\n",nx);
130 // printf(" ny [number of Y combination]:.. %i\n",ny);
131 printf("========================================================================\n");
132 //
133 };
134
135 void CaloLat::Draw(){
136 //
137 Process();
138 Draw(-1,-1);
139 };
140
141 void Calo2D::Draw(){
142 //
143 Process();
144 Draw(-1);
145 };
146
147 void CaloLat::Draw(Int_t view,Int_t plane){
148 //
149 Int_t minv = 0;
150 Int_t maxv = 0;
151 Int_t minp = 0;
152 Int_t maxp = 0;
153 //
154 if ( view == -1 ){
155 minv = 0;
156 maxv = 2;
157 } else {
158 minv = view;
159 maxv = view+1;
160 };
161
162 if ( plane == -1 ){
163 minp = 0;
164 maxp = 22;
165 } else {
166 minp = plane;
167 maxp = plane+1;
168 };
169 //
170 Process();
171 //
172 gStyle->SetLabelSize(0.04);
173 gStyle->SetNdivisions(510,"XY");
174 //
175 for (Int_t v=minv; v<maxv;v++){
176 for (Int_t p=minp; p<maxp;p++){
177 TString hid = Form("clatv%ip%i",v,p);
178 TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
179 if ( tc ){
180 // tc->Clear();
181 } else {
182 tc = new TCanvas(hid,hid);
183 };
184 //
185 TString thid = Form("hlatv%ip%i",v,p);
186 TH1F *th = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
187 if ( th ) th->Delete();
188 // th->Clear();
189 // th->Reset();
190 // } else {
191 th = new TH1F(thid,thid,96,-0.5,95.5);
192 // };
193 tc->cd();
194 //
195 for (Int_t st=0;st<96;st++){
196 th->Fill(st,estrip[v][p][st]);
197 };
198 th->Draw();
199 tc->Modified();
200 tc->Update();
201 };
202 };
203 //
204 gStyle->SetLabelSize(0);
205 gStyle->SetNdivisions(1,"XY");
206 //
207 };
208
209 void Calo2D::Draw(Int_t plane){
210 //
211 Int_t minp = 0;
212 Int_t maxp = 0;
213 //
214 if ( plane == -1 ){
215 minp = 0;
216 maxp = 23;
217 } else {
218 minp = plane;
219 maxp = plane+1;
220 };
221 //
222 Process();
223 //
224 gStyle->SetLabelSize(0.04);
225 gStyle->SetNdivisions(510,"XY");
226 //
227 for (Int_t p=minp; p<maxp;p++){
228 TString hid = Form("c2dp%i",p);
229 TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
230 if ( tc ){
231 // tc->Clear();
232 } else {
233 tc = new TCanvas(hid,hid);
234 };
235 //
236 TString thid = Form("h2dp%i",p);
237 TH2F *th = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
238 if ( th ) th->Delete();
239 // th->Clear();
240 // th->Reset();
241 // } else {
242 Int_t minx = smax[p] - 10;
243 if ( minx < 0 ) minx = 0;
244 Int_t maxx = minx + 20;
245 if ( maxx > 95 ){
246 maxx = 95;
247 minx = 75;
248 };
249 Int_t miny = smay[p] - 10;
250 if ( miny < 0 ) miny = 0;
251 Int_t maxy = miny + 20;
252 if ( maxy > 95 ){
253 maxy = 95;
254 miny = 75;
255 };
256 th = new TH2F(thid,thid,20,(Float_t)minx-0.5,(Float_t)maxx-0.5,20,(Float_t)miny-0.5,(Float_t)maxy-0.5);
257 // th = new TH2F(thid,thid,96,-0.5,95.5,96,-0.5,95.5);
258 // };
259 tc->cd();
260 //
261 for (Int_t stx=minx;stx<maxx+1;stx++){
262 for (Int_t sty=miny;sty<maxy+1;sty++){
263 th->Fill(stx,sty,estrip[p][stx][sty]);
264 };
265 };
266 gStyle->SetPalette(1);
267 // tc->SetLogz();
268 // th->Draw("colbox");
269 th->Draw("cont4");
270 tc->Modified();
271 tc->Update();
272 };
273 //
274 gStyle->SetLabelSize(0);
275 gStyle->SetNdivisions(1,"XY");
276 //
277 };
278
279 void CaloLat::Delete(){
280 Clear();
281 //delete this;
282 };
283
284 void Calo2D::Delete(){
285 Clear();
286 //delete this;
287 };
288
289
290 void CaloLat::Process(){
291 //
292 if ( !L2 ){
293 printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
294 printf(" ERROR: CaloHough variables not filled \n");
295 return;
296 };
297 //
298 Bool_t newentry = false;
299 //
300 if ( L2->IsORB() ){
301 if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
302 newentry = true;
303 OBT = L2->GetOrbitalInfo()->OBT;
304 PKT = L2->GetOrbitalInfo()->pkt_num;
305 atime = L2->GetOrbitalInfo()->absTime;
306 };
307 } else {
308 newentry = true;
309 };
310 //
311 if ( !newentry ) return;
312 //
313 if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
314 //
315 Clear();
316 //
317 // let's start
318 //
319 memset(estrip,0, 4224*sizeof(Float_t));
320 Float_t mip1 = 0.;
321 Int_t view1 = 0;
322 Int_t plane1 = 0;
323 Int_t strip1 = 0;
324 //
325 for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
326 mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
327 estrip[view1][plane1][strip1] = mip1;
328 };
329 //
330 if ( debug ) this->Print();
331 if ( debug ) printf(" exit \n");
332 //
333 };
334
335 void Calo2D::Process(){
336 //
337 if ( !L2 ){
338 printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
339 printf(" ERROR: CaloHough variables not filled \n");
340 return;
341 };
342 //
343 Bool_t newentry = false;
344 //
345 if ( L2->IsORB() ){
346 if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
347 newentry = true;
348 OBT = L2->GetOrbitalInfo()->OBT;
349 PKT = L2->GetOrbitalInfo()->pkt_num;
350 atime = L2->GetOrbitalInfo()->absTime;
351 };
352 } else {
353 newentry = true;
354 };
355 //
356 if ( !newentry ) return;
357 //
358 if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
359 //
360 Clear();
361 //
362 // let's start
363 //
364 Float_t es[2][22][96];
365 memset(es,0, 4224*sizeof(Float_t));
366 memset(estrip,0, 4224*sizeof(Float_t));
367 Float_t mip1 = 0.;
368 Int_t view1 = 0;
369 Int_t plane1 = 0;
370 Int_t strip1 = 0;
371 //
372 for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
373 mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
374 es[view1][plane1][strip1] = mip1;
375 };
376 //
377 Int_t plane2 = 0;
378 Float_t emax[23];
379 memset(emax,0,sizeof(Float_t)*23);
380 memset(smax,0,sizeof(Int_t)*23);
381 memset(smay,0,sizeof(Int_t)*23);
382 //
383 //
384 for (Int_t p=0; p < 23 ; p++){
385 //
386 plane1 = p-1;
387 plane2 = p;
388 //
389 if ( p == 0 ){
390 plane1 = -1;
391 plane2 = 0;
392 };
393 if ( p == 22 ){
394 plane1 = 21;
395 plane2 = -1;
396 };
397 //
398 for (Int_t s=0; s < 96 ; s++){ // x
399 for (Int_t ss=0; ss < 96 ; ss++){ // y
400 if ( p == 0 ){
401 estrip[p][s][ss] += es[1][plane2][ss];
402 if ( (es[1][plane2][ss]) > emax[p] ){
403 smax[p] = 45;
404 smay[p] = ss;
405 emax[p] = es[1][plane2][ss] ;
406 };
407 };
408 if ( p > 0 && p < 22 ){
409 estrip[p][s][ss] += es[0][plane1][s] + es[1][plane2][ss];
410 if ( (es[0][plane1][s] + es[1][plane2][ss]) > emax[p] ){
411 smax[p] = s;
412 smay[p] = ss;
413 emax[p] = es[0][plane1][s] + es[1][plane2][ss] ;
414 };
415 };
416 if ( p == 22 ){
417 estrip[p][s][ss] += es[0][plane1][s];
418 if ( (es[1][plane2][s]) > emax[p] ){
419 smax[p] = s;
420 smay[p] = 45;
421 emax[p] = es[1][plane2][s] ;
422 };
423 };
424 };
425 };
426 //
427 };
428 //
429 if ( debug ) this->Print();
430 if ( debug ) printf(" exit \n");
431 //
432 };
433
434
435 /**
436 * Default constructor
437 */
438 CaloLong::CaloLong(){
439 Clear();
440 };
441
442 CaloLong::CaloLong(PamLevel2 *l2p){
443 //
444 Clear();
445 //
446 L2 = l2p;
447 //
448 if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
449 //
450 OBT = 0;
451 PKT = 0;
452 atime = 0;
453 //
454 clp = L2->GetCaloLevel2();
455 //
456 sel = true;
457 cont = false;
458 N = 0;
459 NC = 22;
460 mask18b = -1;
461 //
462 no18x = true;
463 debug = false;
464 maskXE = false;
465 maskXO = false;
466 maskYE = false;
467 maskYO = false;
468 //
469 };
470
471 void CaloLong::MaskSection(TString sec){
472 sec.ToUpper();
473 if ( sec.Contains("XO") ) maskXO = true;
474 if ( sec.Contains("YO") ) maskYO = true;
475 if ( sec.Contains("XE") ) maskXE = true;
476 if ( sec.Contains("YE") ) maskYE = true;
477 }
478
479 void CaloLong::UnMaskSections(){
480 this->UnMaskSection("XEXOYEYO");
481 }
482
483 void CaloLong::UnMaskSection(TString sec){
484 sec.ToUpper();
485 if ( sec.Contains("XO") ) maskXO = false;
486 if ( sec.Contains("YO") ) maskYO = false;
487 if ( sec.Contains("XE") ) maskXE = false;
488 if ( sec.Contains("YE") ) maskYE = false;
489 }
490
491 void CaloLong::Clear(){
492 //
493 memset(eplane,0, 2*22*sizeof(Float_t));
494 //
495 chi2 = 0.;
496 ndf = 0.;
497 E0 = 0.;
498 a = 0.;
499 b = 0.;
500 errE0 = 0.;
501 erra = 0.;
502 errb = 0.;
503 etmax = 0.;
504 asymm = 0.;
505 fitresult = 0;
506 //
507 X0pl = 0.76;
508 //
509 };
510
511 void CaloLong::Print(){
512 //
513 Fit();
514 //
515 printf("==================== Calorimeter Longitudinal Profile =======================\n");
516 printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
517 printf(" fitresult:.. %i\n",fitresult);
518 printf(" chi2 :.. %f\n",chi2);
519 printf(" ndf :.. %f\n",ndf);
520 printf(" nchi2 :.. %f\n",chi2/ndf);
521 printf(" E0 :.. %f\n",E0);
522 printf(" E0/260. :.. %f\n",E0/260.);
523 printf(" a :.. %f\n",a);
524 printf(" b :.. %f\n",b);
525 printf(" errE0 :.. %f\n",errE0);
526 printf(" erra :.. %f\n",erra);
527 printf(" errb :.. %f\n",errb);
528 printf(" asymm :.. %f\n",asymm);
529 printf(" tmax :.. %f\n",((a-1.)/b));
530 printf(" etmax :.. %f\n",etmax);
531 printf(" X0pl :.. %f\n",X0pl);
532 printf("========================================================================\n");
533 //
534 };
535
536 void CaloLong::SetNoWpreSampler(Int_t n){
537 //
538 if ( NC+n <= 22 && NC+n >= 0 ){
539 N = n;
540 } else {
541 printf(" ERROR! Calorimeter is made of 22 W planes\n");
542 printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
543 printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
544 NC = 22;
545 N = 0;
546 };
547 }
548
549 void CaloLong::SetNoWcalo(Int_t n){
550 if ( N+n <= 22 && N+n >= 0 ){
551 NC = n;
552 } else {
553 printf(" ERROR! Calorimeter is made of 22 W planes\n");
554 printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
555 printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
556 NC = 22;
557 N = 0;
558 };
559 }
560
561 void CaloLong::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){
562 this->SetNoWpreSampler(0);
563 this->SetNoWcalo(0);
564 if ( NoWpreSampler < NoWcalo ){
565 this->SetNoWpreSampler(NoWpreSampler);
566 this->SetNoWcalo(NoWcalo);
567 } else {
568 this->SetNoWcalo(NoWcalo);
569 this->SetNoWpreSampler(NoWpreSampler);
570 };
571 }
572
573 void CaloLong::Process(){
574 //
575 if ( !L2 ){
576 printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
577 printf(" ERROR: CaloHough variables not filled \n");
578 return;
579 };
580 //
581 Bool_t newentry = false;
582 //
583 if ( L2->IsORB() ){
584 if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
585 newentry = true;
586 OBT = L2->GetOrbitalInfo()->OBT;
587 PKT = L2->GetOrbitalInfo()->pkt_num;
588 atime = L2->GetOrbitalInfo()->absTime;
589 };
590 } else {
591 newentry = true;
592 };
593 //
594 if ( !newentry ) return;
595 //
596 if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
597 //
598 Clear();
599 //
600 // let's start
601 //
602 if ( cont ){
603 for (Int_t i=0; i<22; i++){
604 if ( i == (18+N) ){
605 mask18b = 18 + N;
606 break;
607 };
608 };
609 };
610 //
611 if ( sel ){
612 for (Int_t i=0; i<22; i++){
613 if ( i == (18-N) ){
614 mask18b = 18 - N;
615 break;
616 };
617 };
618 };
619 //
620 // if ( mask18b == 18 ) mask18b = -1;
621 //
622 Int_t view = 0;
623 Int_t plane = 0;
624 Int_t strip = 0;
625 Float_t mip = 0.;
626 Bool_t gof = true;
627 for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){
628 mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
629 gof = true;
630 if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
631 if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
632 if ( maskYE && (plane%2)!=0 && view==0 ) gof = false;
633 if ( maskYO && (plane%2)==0 && view==0 ) gof = false;
634 if ( gof ) eplane[view][plane] += mip;
635 };
636 //
637 // inclination factor (stolen from Daniele's code)
638 //
639 Float_t ytgx = 0;
640 Float_t ytgy = 0;
641 ytgx = 0.76 * clp->tanx[0];
642 ytgy = 0.76 * clp->tany[0];
643 X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
644 //
645 // Find experimental plane of maximum
646 //
647 Int_t pmax = 0;
648 Int_t vmax = 0;
649 Float_t emax = 0.;
650 for (Int_t v=0; v<2; v++){
651 for (Int_t i=0; i<22; i++){
652 if ( eplane[v][i] > emax ){
653 emax = eplane[v][i];
654 vmax = v;
655 pmax = i;
656 };
657 };
658 };
659 //
660 //
661 //
662 if ( vmax == 0 ) pmax++;
663 etmax = pmax * X0pl;
664 //
665 if ( debug ) this->Print();
666 if ( debug ) printf(" exit \n");
667 //
668 };
669
670 void CaloLong::SetEnergies(Float_t myene[][22]){
671 //
672 if ( !L2 ){
673 printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
674 printf(" ERROR: CaloHough variables not filled \n");
675 return;
676 };
677 //
678 Bool_t newentry = false;
679 //
680 if ( L2->IsORB() ){
681 if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
682 newentry = true;
683 OBT = L2->GetOrbitalInfo()->OBT;
684 PKT = L2->GetOrbitalInfo()->pkt_num;
685 atime = L2->GetOrbitalInfo()->absTime;
686 };
687 } else {
688 newentry = true;
689 };
690 //
691 if ( !newentry ) return;
692 //
693 if ( debug ) printf(" SET ENERGIES Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
694 //
695 Clear();
696 //
697 // let's start
698 //
699 if ( cont ){
700 for (Int_t i=0; i<22; i++){
701 if ( i == (18+N) ){
702 mask18b = 18 + N;
703 break;
704 };
705 };
706 };
707 //
708 if ( sel ){
709 for (Int_t i=0; i<22; i++){
710 if ( i == (18-N) ){
711 mask18b = 18 - N;
712 break;
713 };
714 };
715 };
716 //
717 // if ( mask18b == 18 ) mask18b = -1;
718 //
719 Int_t view = 0;
720 Int_t plane = 0;
721 Bool_t gof = true;
722 for (view=0; view < 2; view++){
723 for (plane=0; plane < 22; plane++){
724 gof = true;
725 if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
726 if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
727 if ( maskYE && (plane%2)!=0 && view==0 ) gof = false;
728 if ( maskYO && (plane%2)==0 && view==0 ) gof = false;
729 if ( gof ) eplane[view][plane] = myene[view][plane];
730 if ( debug ) printf(" XE %i XO %i YE %i YO %i eplane %i %i = %f ........... myene %i %i = %f\n", maskXE, maskXO, maskYE, maskYO,view,plane,eplane[view][plane],view,plane,myene[view][plane]);
731 };
732 };
733 //
734 // inclination factor (stolen from Daniele's code)
735 //
736 Float_t ytgx = 0;
737 Float_t ytgy = 0;
738 ytgx = 0.76 * clp->tanx[0];
739 ytgy = 0.76 * clp->tany[0];
740 X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
741 //
742 // Find experimental plane of maximum
743 //
744 Int_t pmax = 0;
745 Int_t vmax = 0;
746 Float_t emax = 0.;
747 for (Int_t v=0; v<2; v++){
748 for (Int_t i=0; i<22; i++){
749 if ( eplane[v][i] > emax ){
750 emax = eplane[v][i];
751 vmax = v;
752 pmax = i;
753 };
754 };
755 };
756 //
757 //
758 //
759 if ( vmax == 0 ) pmax++;
760 etmax = pmax * X0pl;
761 //
762 if ( debug ) this->Print();
763 if ( debug ) printf(" exit \n");
764 //
765 };
766
767 Double_t ccurve(Double_t *ti,Double_t *par){
768 //
769 Double_t t = *ti;
770 Double_t cE0 = par[0];
771 Double_t ca = par[1];
772 Double_t cb = par[2];
773 Double_t gammaa = TMath::Gamma(ca);
774 //
775 Double_t value = cE0 * cb * ( pow((cb * t),(ca - 1)) * exp( -cb * t ) ) / gammaa;
776 //
777 return value;
778 //
779 }
780
781 void CaloLong::Fit(){
782 this->Fit(false);
783 };
784
785 void CaloLong::Fit(Bool_t draw){
786 //
787 Process();
788 //
789 //
790 if ( !L2 ){
791 printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
792 printf(" ERROR: CaloHough variables not filled \n");
793 return;
794 };
795 //
796 Bool_t newentry = false;
797 //
798 if ( L2->IsORB() ){
799 if ( L2->GetOrbitalInfo()->pkt_num != fPKT || L2->GetOrbitalInfo()->OBT != fOBT || L2->GetOrbitalInfo()->absTime != fatime ){
800 newentry = true;
801 fOBT = L2->GetOrbitalInfo()->OBT;
802 fPKT = L2->GetOrbitalInfo()->pkt_num;
803 fatime = L2->GetOrbitalInfo()->absTime;
804 };
805 } else {
806 newentry = true;
807 };
808 //
809 if ( !newentry ) return;
810 //
811 if ( debug ) printf(" Start fitting event at OBT %u PKT %u time %u \n",fOBT,fPKT,fatime);
812 //
813 if ( draw ){
814 gStyle->SetLabelSize(0.04);
815 gStyle->SetNdivisions(510,"XY");
816 };
817 //
818 TString hid = Form("clongfit");
819 TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
820 // if ( tc ) tc->Delete();
821 // if ( tc ) tc->Close();
822 if ( !tc && draw ){
823 tc = new TCanvas(hid,hid);
824 } else {
825 if ( tc ) tc->cd();
826 };
827 //
828 TString thid = Form("hlongfit");
829 TH1F *th = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
830 if ( th ) th->Delete();
831 Float_t xpos = 0.;
832 Float_t enemip = 0.;
833 Float_t xmax = NC * X0pl + 0.2;
834 // th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax);
835 th = new TH1F(thid,thid,100,-0.2,xmax);
836 //
837 // AGH, BUG!
838 //
839 Int_t mmin = 0;
840 Int_t mmax = 0;
841 if ( cont ){
842 mmin = N;
843 mmax = NC+N;
844 } else {
845 mmin = 0;
846 mmax = NC;
847 };
848 //
849 Float_t qtotparz = 0.;
850 for (Int_t st=mmin;st<mmax+1;st++){
851 enemip = 0.;
852 xpos = (st - mmin) * X0pl;
853 if ( st > mmin && st < mmax ){
854 if ( no18x && ( st == 18+1 || st == mask18b+1 )){
855 if ( !maskYO ){
856 enemip = 2. * eplane[1][st];
857 } else {
858 enemip = eplane[1][st];
859 };
860 } else {
861 enemip = eplane[0][st-1] + eplane[1][st];
862 };
863 } else {
864 if ( st == mmin ){
865 if ( !maskYE ){
866 enemip = 2. * eplane[1][st];
867 } else {
868 enemip = eplane[1][st];
869 };
870 };
871 if ( st == mmax ){
872 if ( !maskXE ){
873 enemip = 2. * eplane[0][st-1];
874 } else {
875 enemip = eplane[0][st-1];
876 };
877 };
878 };
879 //
880 qtotparz += enemip;
881 if ( enemip > 0. ){
882 th->Fill(xpos,enemip);
883 if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
884 };
885 //
886 // for (Int_t v=1; v>=0;v--)// {
887 // //
888 // if ( v == 1 ){
889 // xpos = (st - N) * X0pl;
890 // } else {
891 // xpos = (st + 1 - N) * X0pl;
892 // };
893 // //
894 // if ( no18x && st == 18 && v == 0 ){
895 // // skip plane 18x
896 // } else {
897 // if ( v == 1 && st == mask18b ){
898 // // emulate plane 18x
899 // } else {
900 // if ( eplane[v][st] > 0. ){
901 // th->Fill(xpos,eplane[v][st]);
902 // if ( debug ) printf(" Filling: st %i v %i xpos %f energy %f \n",st,v,xpos,eplane[v][st]);
903 // };
904 // };
905 // };
906 // //
907 // };
908 };
909 //
910 TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);
911 if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz);
912 E0 = qtotparz;
913 // E0 = clp->qtot;
914 a = 5.;
915 b = 0.5;
916 if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b);
917 lfit->SetParameters(E0,a,b);
918 // lfit->SetParLimits(0,0.,1000.);
919 // lfit->SetParLimits(1,-1.,80.);
920 // lfit->SetParLimits(2,-1.,10.);
921 TString optio;
922 if ( debug ){
923 // optio = "MERBOV";
924 // optio = "MEROV";
925 // optio = "EROV";
926 optio = "RNOV";
927 if ( draw ) optio = "ROV";
928 } else {
929 // optio = "MERNOQ";
930 // optio = "ERNOQ";
931 optio = "RNOQ";
932 if ( draw ) optio = "ROQ";
933 };
934 //
935 if ( debug ) printf(" OK, start the fitting procedure...\n");
936 //
937 fitresult = th->Fit("lfit",optio);
938 //
939 if ( debug ) printf(" the fit is done! result: %i \n",fitresult);
940 //
941 E0 = lfit->GetParameter(0);
942 a = lfit->GetParameter(1);
943 b = lfit->GetParameter(2);
944 errE0 = lfit->GetParError(0);
945 erra = lfit->GetParError(1);
946 errb = lfit->GetParError(2);
947 chi2 = lfit->GetChisquare();
948 ndf = lfit->GetNDF();
949 Float_t tmax = 0.;
950 if ( debug ) printf(" Parameters are retrieved \n");
951 if ( b != 0 ) tmax = (a - 1.)/b;
952 //
953 if ( fitresult != 0 ){
954 if ( debug ) printf(" The fit failed, no integrals calculation and asymm is set to -1. \n");
955 asymm = -1.;
956 } else {
957 Int_t npp = 1000;
958 double *xp=new double[npp];
959 double *wp=new double[npp];
960 lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12);
961 Float_t imax = lfit->IntegralFast(npp,xp,wp,0.,tmax);
962 // Float_t imax = lfit->Integral(0.,tmax);
963 if ( debug ) printf(" Integral till maximum (%f): %f \n",tmax,imax);
964 Int_t np = 1000;
965 double *x=new double[np];
966 double *w=new double[np];
967 lfit->CalcGaussLegendreSamplingPoints(np,x,w,1e-12);
968 Float_t i10max = lfit->IntegralFast(np,x,w,0.,10.*tmax);
969 delete x;
970 delete w;
971 delete xp;
972 delete wp;
973 // Float_t i10max = lfit->Integral(0.,10.*tmax);
974 if ( debug ) printf(" Integral: %f \n",i10max);
975 //
976 if ( i10max != imax ){
977 asymm = imax / (i10max-imax);
978 } else {
979 if ( debug ) printf(" i10max == imax, asymm undefined\n");
980 asymm = -2.;
981 };
982 if ( asymm != asymm ){
983 if ( debug ) printf(" asymm is nan \n");
984 asymm = -3.;
985 };
986 //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax));
987 if ( debug ) printf(" Asymmetry has been calculated \n");
988 };
989 //
990 if ( asymm < 0. || ndf <= 0. || chi2 < 0. || tmax < 0. ){
991 if ( debug ) printf(" Funny asymm||ndf||chi2||tmax values, fit failed \n");
992 fitresult = 100;
993 };
994 //
995 if ( draw ){
996 //
997 tc->cd();
998 // gStyle->SetOptStat(11111);
999 tc->SetTitle();
1000 th->SetTitle("");
1001 th->SetName("");
1002 th->SetMarkerStyle(20);
1003 // axis titles
1004 th->SetXTitle("Depth [X0]");
1005 th->SetYTitle("Energy [MIP]");
1006 th->DrawCopy("Perror");
1007 lfit->Draw("same");
1008 tc->Modified();
1009 tc->Update();
1010 //
1011 gStyle->SetLabelSize(0);
1012 gStyle->SetNdivisions(1,"XY");
1013 //
1014 } else {
1015 if ( th ) th->Delete();
1016 };
1017 //
1018 delete lfit;
1019 //
1020 };
1021
1022 void CaloLong::Draw(){
1023 //
1024 Process();
1025 Draw(-1);
1026 };
1027
1028 void CaloLong::Draw(Int_t view){
1029 //
1030 Int_t minv = 0;
1031 Int_t maxv = 0;
1032 //
1033 if ( view == -1 ){
1034 maxv = -1;
1035 } else {
1036 minv = view;
1037 maxv = view+1;
1038 };
1039 //
1040 Process();
1041 //
1042 gStyle->SetLabelSize(0.04);
1043 gStyle->SetNdivisions(510,"XY");
1044 //
1045 if ( maxv != -1 ){
1046 for (Int_t v=minv; v<maxv;v++){
1047 TString hid = Form("clongv%i",v);
1048 TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1049 if ( tc ){
1050 // tc->Clear();
1051 } else {
1052 tc = new TCanvas(hid,hid);
1053 };
1054 //
1055 TString thid = Form("hlongv%i",v);
1056 TH1F *th = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1057 if ( th ) th->Delete();
1058 // th->Clear();
1059 // th->Reset();
1060 // } else {
1061 th = new TH1F(thid,thid,22,-0.5,21.5);
1062 // };
1063 tc->cd();
1064 //
1065 for (Int_t st=0;st<22;st++){
1066 th->Fill(st,eplane[v][st]);
1067 };
1068 th->Draw();
1069 tc->Modified();
1070 tc->Update();
1071 };
1072 } else {
1073 //
1074 TString hid = Form("clongvyvx");
1075 TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1076 if ( tc ){
1077 } else {
1078 tc = new TCanvas(hid,hid);
1079 };
1080 //
1081 TString thid = Form("hlongvyvx");
1082 TH1F *th = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1083 if ( th ) th->Delete();
1084 th = new TH1F(thid,thid,44,-0.5,43.5);
1085 tc->cd();
1086 Int_t pp=0;
1087 for (Int_t st=0;st<22;st++){
1088 for (Int_t v=1; v>=0;v--){
1089 //
1090 th->Fill(pp,eplane[v][st]);
1091 //
1092 pp++;
1093 };
1094 };
1095 th->Draw();
1096 tc->Modified();
1097 tc->Update();
1098 };
1099 //
1100 gStyle->SetLabelSize(0);
1101 gStyle->SetNdivisions(1,"XY");
1102 //
1103 };
1104
1105 void CaloLong::Delete(){
1106 Clear();
1107 //delete this;
1108 };

  ViewVC Help
Powered by ViewVC 1.1.23