/[PAMELA software]/PamCut/CaloAxis2.h
ViewVC logotype

Contents of /PamCut/CaloAxis2.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Wed May 27 13:30:01 2009 UTC (15 years, 6 months ago) by pam-fi
Branch: DEV
CVS Tags: v0r00
Changes since 1.1: +0 -0 lines
File MIME type: text/plain
First import

1 #ifndef CALOAXIS2_H_
2 #define CALOAXIS2_H_
3
4 #if !defined(__CINT__) || defined(__MAKECINT__)
5
6
7 #include <TString.h>
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TMath.h>
11 #include <TLine.h>
12 #include <TPolyMarker.h>
13 #include <TSelector.h>
14 #include <TFile.h>
15
16 #include <stdlib.h>
17 #include <iostream>
18 using namespace std;
19
20 #include <PamLevel2.h>
21 //#include <CaloEvent.h>
22
23 #endif // !defined(__CINT__) || defined(__MAKECINT__)
24 //===============================================================================
25 //
26 //
27 //
28 //
29 //===============================================================================
30 class CaloEvent{
31
32 public:
33 // calo event
34 // ------------------------------
35 // hit strips
36 // ------------------------------
37 vector<int> ids; //strip 0-95
38 vector<int> idp; //plane 0-21
39 vector<float> q; // signal
40 vector<float> zcoord;
41 vector<float> xycoord;
42
43
44 // ------------------------------
45 // methods
46 // ------------------------------
47 void Add(int iis, int iip, float fq, float fz, float fxy){
48 ids.push_back(iis);
49 idp.push_back(iip);
50 q.push_back(fq);
51 zcoord.push_back(fz);
52 xycoord.push_back(fxy);
53 };
54
55
56
57 int GetN(){return (int)(q.size());};
58
59 void Reset(){
60 ids.clear();
61 idp.clear();
62 q.clear();
63 zcoord.clear();
64 xycoord.clear();
65 };
66
67 void Delete(){ Reset(); };
68
69 ~CaloEvent(){ Delete(); };
70
71 CaloEvent(){ };
72
73 CaloEvent( CaloLevel1* calo, int view, Bool_t usemechanicalalignement){ Set(calo,view,usemechanicalalignement); };
74
75 CaloEvent( CaloLevel1* calo, int view){ Set(calo,view); };
76
77 void Set(CaloLevel1* calo, int view, Bool_t usemechanicalalignement){
78
79 Reset();
80 // cout << " CaloEvent::Set()"<<endl;
81 CaloStrip cstrip;
82 cstrip = CaloStrip(calo,usemechanicalalignement);
83
84 for(int ih=0;ih<calo->istrip;ih++){
85 int iv=-1;
86 int ip=-1;
87 int is=-1;
88 calo->DecodeEstrip(ih,iv,ip,is);
89 if(iv==view){
90 cstrip.Set(view,ip,is);
91 float fxy = 0.;
92 if( !view ) fxy = cstrip.GetX();
93 else fxy = cstrip.GetY();
94 float fz = cstrip.GetZ();
95 float fq = cstrip.GetE();
96 if(fq>0)Add(is,ip,fq,fz,fxy);
97 }
98 }
99
100 };
101 void Set(CaloLevel1* calo, int view){ Set(calo,view,0); };
102
103
104 };
105 //===============================================================================
106 //
107 //
108 //
109 //
110 //===============================================================================
111 class CaloAxis {
112
113 public:
114
115 CaloEvent *cevent;
116
117 // ------------------------------
118 // selected points along the axis
119 // ------------------------------
120 // fitted points
121 vector<float> x;// independent coordinates (z)
122 vector<float> y;// dependente coordinate (x/y)
123 vector<float> w;// fit weight
124 vector<float> q;// signal
125
126 // ------------------------------
127 // axis parameters
128 // ------------------------------
129 Float_t par[2];
130 Float_t covpar[2][2];
131
132 // ------------------------------
133 // general variables
134 // ------------------------------
135 float cog[22];// baricenter
136 float qpl[22]; // charge in each plane
137
138 // ------------------------------
139 // objects to draw the axis
140 // ------------------------------
141 TPolyMarker *sel;
142 TLine *axis;
143
144
145 // ------------------------------
146 // methods
147 // ------------------------------
148 void Init(){
149 for(Int_t i=0; i<2;i++){
150 par[i]=0;
151 for(Int_t j=0; j<2;j++)covpar[i][j]=0;
152 };
153 sel = 0;
154 axis = 0;
155 for(Int_t i=0; i<22; i++){
156 cog[i]=0;
157 qpl[i]=0;
158 }
159 // if(cevent)cevent->Reset();
160 }
161
162 void Set(CaloLevel1* calo, Int_t view){ cevent->Set(calo,view,0); };
163 void Set(CaloLevel1* calo, Int_t view, Bool_t usemechanicalalignement){ cevent->Set(calo,view,usemechanicalalignement); };
164
165 CaloAxis(){
166 cevent = new CaloEvent();
167 Init();
168 };
169
170 CaloAxis(CaloLevel1* calo, Int_t view){
171 cevent = new CaloEvent();
172 Init();
173 Set(calo,view);
174 };
175 CaloAxis(CaloLevel1* calo, Int_t view, Bool_t usemechanicalalignement){
176 cevent = new CaloEvent();
177 Init();
178 Set(calo,view,usemechanicalalignement);
179 };
180
181
182 void Reset(){
183 // cout << " CaloAxis::Reset()"<<endl;
184 x.clear();
185 y.clear();
186 w.clear();
187 q.clear();
188 // cout << sel << endl;
189 if(sel)sel->Delete();
190 // cout << axis << endl;
191 if(axis)axis->Delete();
192 // if(cevent)cevent->Reset();
193 Init();
194 }
195
196 void Delete(){ Reset(); delete cevent; };
197
198 ~CaloAxis(){ Delete(); };
199
200
201 void Add(float xin, float yin){
202 x.push_back(xin);
203 y.push_back(yin);
204 w.push_back(1.);
205 q.push_back(0.);
206 /* Int_t nin = x.size(); */
207 /* if(nin>2)Fit(); */
208 };
209 void Add(float xin, float yin, float qin){
210 x.push_back(xin);
211 y.push_back(yin);
212 q.push_back(qin);
213 w.push_back(1.);
214 /* Int_t nin = x.size(); */
215 /* if(nin>2)Fit(); */
216 };
217 void Add(float xin, float yin, float qin, float win){
218 x.push_back(xin);
219 y.push_back(yin);
220 w.push_back(win);
221 q.push_back(qin);
222 /* Int_t nin = x.size(); */
223 /* if(nin>2)Fit(); */
224 };
225
226 int Fit(){
227
228 // cout << " CaloAxis::Fit()"<<endl;
229 Float_t SSS = 0;
230 Float_t SSX = 0;
231 Float_t SXX = 0;
232
233 Float_t SSY = 0;
234 Float_t SXY = 0;
235
236 Int_t np=0;
237 for(Int_t i=0; i<(int)(x.size()); i++){
238 SSS += w[i]*w[i];
239 SSX += x[i]*w[i]*w[i];
240 SXX += x[i]*x[i]*w[i]*w[i];
241 SSY += y[i]*w[i]*w[i];
242 SXY += x[i]*y[i]*w[i]*w[i];
243 if(w[i])np++;
244 }
245 Float_t det = SSS*SXX-SSX*SSX;
246 if(det==0)return 0;
247 if(np<3)return 0;
248 // cout << np << " points fitted -- ok "<<endl;
249 par[0] = (SSY*SXX-SXY*SSX)/det;
250 par[1] = (SXY*SSS-SSY*SSX)/det;
251 return 1;
252 };
253
254 void Print(){
255 // useful for debug for(Int_t i=0; i<(int)(x.size()); i++)if(w[i])cout << x[i] << " - "<<y[i]<<endl;
256 }
257
258 float GetPar(Int_t i){return par[i];};
259
260 int GetN(){return (int)(x.size());};
261
262 float* GetX(){
263 float *xout = new float[(int)(x.size())];
264 for(Int_t i=0; i<(int)(x.size()); i++)xout[i]=x[i];
265 return xout;
266 }
267
268 float* GetY(){
269 float *yout = new float[(int)(y.size())];
270 for(Int_t i=0; i<(int)(y.size()); i++) yout[i]=y[i];
271 return yout;
272 }
273
274 float* GetQ(){
275 float *qout = new float[(int)(q.size())];
276 for(int i=0; i<(int)(q.size()); i++) qout[i]=q[i];
277 return qout;
278 }
279
280 float GetChi2(){
281 float chi2=0;
282 int nchi=0;
283 for(Int_t i=0; i<(int)(y.size()); i++){
284 chi2 += w[i]*w[i]*(y[i]-par[0]-par[1]*x[i])*(y[i]-par[0]-par[1]*x[i]);
285 nchi += (int)w[i];
286 };
287
288 if(nchi<3) return -999; // original
289 // if(nchi<3 || chi2==0)return -999; // modified by Sergio
290
291 chi2=chi2/(nchi-2);
292 return chi2;
293 }
294
295 float GetQaxis(){
296 float qaxis=0;
297 for(Int_t i=0; i<(int)(q.size()); i++) qaxis +=q[i];
298 // for(Int_t i=0; i<(int)(q.size()); i++) cout<<q[i]<<endl;
299 return qaxis;
300 }
301
302 float GetQaxis(float toll){
303 float qaxis=0;
304 for(int ih=0;ih<cevent->GetN(); ih++){
305 float x = cevent->xycoord[ih];
306 float z = cevent->zcoord[ih];
307 float q = cevent->q[ih];
308 //int ip = cevent->idp[ih];
309 float d = fabs(x-par[0]-par[1]*z)/sqrt(par[1]*par[1]+1);
310 if( d < toll )qaxis+=q;
311 }
312 return qaxis;
313 }
314
315 float GetCOG(Int_t ip){
316 if(ip>=0 && ip<22)return cog[ip];
317 else return 0;
318 }
319
320 float GetQ(Int_t ip){
321 if(ip>=0 && ip<22)return qpl[ip];
322 else return 0;
323 }
324
325 float GetQstrip(){
326
327 float qaxis=0;
328 for(Int_t i=0; i<(int)(q.size()); i++) qaxis +=q[i];
329 if(q.size()==0) return 0.;
330 return qaxis/q.size();
331 }
332
333 float GetYfit(float x){
334 return par[0]+x*par[1];
335 }
336
337 float GetXmin(){
338 float zmin=9999;
339 for(Int_t i=0; i<(int)(x.size()); i++) if(x[i]<zmin) zmin = x[i];
340 return zmin;
341 }
342 float GetXmax(){
343 float zmax=-9999;
344 for(Int_t i=0; i<(int)(x.size()); i++) if(x[i]>zmax) zmax = x[i];
345 return zmax;
346
347 }
348
349
350 void Draw(){
351
352 // cout << " CaloAxis::Draw()"<<endl;
353 if(GetN()<3)return;
354
355 sel = new TPolyMarker(GetN(),GetY(),GetX());
356 sel->SetMarkerSize(0.3);
357 sel->SetMarkerColor(5);
358
359 axis = new TLine( par[0]+GetXmin()*par[1], GetXmin(), par[0]+GetXmax()*par[1], GetXmax() );
360 axis->SetLineWidth(1);
361 axis->SetLineColor(5);
362
363 // cout << sel << endl;
364 // cout << axis << endl;
365
366 sel->Draw();
367 axis->Draw();
368 };
369
370
371 //-----------------------------------------------------------------
372 // fit the axis (optimized for non-interacting patterns)
373 //-----------------------------------------------------------------
374 int FitAxis( CaloLevel1* calo , Int_t view , Float_t toll , Bool_t usemechanicalalignement ){
375
376 // cout << "CaloAxis::FitAxis(...)"<<endl;
377
378 Set(calo,view,usemechanicalalignement);
379
380 // ------------------------------
381 // first : all hits included, w=1
382 // ------------------------------
383 for(int ih=0;ih<cevent->GetN(); ih++){
384 float x = cevent->xycoord[ih];
385 float z = cevent->zcoord[ih];
386 float q = cevent->q[ih];
387 Add(z,x,q);
388 }
389
390 if(GetN()<3)return 0;
391 Fit();
392 // useful for debug cout << "(1)"<<endl;
393 // cout << " n. hits :"<<GetN()<<endl;
394 // cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl;
395 // cout << " Chi2 :"<<GetChi2()<<endl;
396 // cout << " Q axis :"<<GetQaxis()<<endl;
397 // cout << " Q strip :"<<GetQstrip()<<endl;
398
399 float P0;
400 float P1;
401 int dof;
402
403 // ---------------------------
404 // second : iteration
405 // ---------------------------
406 int niter=0; // number of iterations
407 int pok[22];
408
409 // prova
410 // pesi lineari
411 float W22=1.;
412 float W1=10.;
413 float b=(W22-W1)/(22.-1.);// (w22 - w1)/(22-1)
414 float a=W1-b;
415
416 for(int i=0; i<22; i++)pok[i]=0;
417 do{
418 niter++;
419 float ttoll = toll; //tolerance (cm)
420 if(niter==1)ttoll = 10.*toll;
421 dof = GetN(); //previous fit
422 P0 = par[0];
423 P1 = par[1];
424 Reset(); //fit reset
425 for(int ih=0;ih<cevent->GetN(); ih++){//loop over selected hits
426 float x = cevent->xycoord[ih];
427 float z = cevent->zcoord[ih];
428 float q = cevent->q[ih];
429 int ip = cevent->idp[ih];
430 float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1);//distance to axis
431 // cout<<d<<endl;
432 if( d < ttoll && (niter==1 || (niter>1 && pok[ip]==1)) ){
433 // if( d < ttoll ){
434 float W=a+b*(ip+1);
435 // cout << ip << " "<<W<<endl;
436 Add(z,x,q,W);
437 pok[ip]=1;
438 }
439 }
440 // break the track if more than 3 planes are missing
441 int nmissing = 0;
442 for(int ip=0; ip<22; ip++){
443 if(pok[ip]==0){
444 nmissing++;
445 }else{
446 nmissing = 0;
447 }
448 if(nmissing==3){
449 for(int ipp=ip+1; ipp<22; ipp++)pok[ipp]=0;
450 break;
451 }
452 }
453
454
455 if(niter==100)break;
456 if(GetN()<3)break;
457 Fit();
458 }while(niter==1 || GetN() != dof);
459 // cout << " >>> "<<GetN()<<endl;
460 if(GetN()<3)return 0;
461 Fit();
462 // cout << "(2)"<<endl;
463 // cout << " n. hits :"<<GetN()<<endl;
464 // cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl;
465 // cout << " Chi2 :"<<GetChi2()<<endl;
466 // cout << " Q axis :"<<GetQaxis()<<endl;
467 // cout << " Q strip :"<<GetQstrip()<<endl;
468
469 // ---------------------------------------------
470 // third : selecting only closest strip-clusters
471 // and fit baricenters
472 // ---------------------------------------------
473 P0 = par[0];
474 P1 = par[1];
475 Reset();
476
477 float dmin[22];
478 int closest[22];
479 for(int ip=0; ip<22; ip++){
480 dmin[ip]=999;
481 closest[ip]=-1;
482 }
483 for(int ih=0;ih<cevent->GetN(); ih++){
484 float x = cevent->xycoord[ih];
485 float z = cevent->zcoord[ih];
486 //float q = cevent->q[ih];
487 int ip = cevent->idp[ih];
488 //int is = cevent->ids[ih];
489 float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1);
490 if( d < toll && d<dmin[ip] && pok[ip]==1 ){
491 // Add(z,x,q,q);
492 closest[ip]=ih;
493 // cog[ip] += x*q;
494 // qpl[ip] += q;
495 }
496 }
497 for(Int_t ip=0; ip<22; ip++){
498 if(closest[ip]!=-1){
499 float x = cevent->xycoord[closest[ip]];
500 float z = cevent->zcoord[closest[ip]];
501 float q = cevent->q[closest[ip]];
502 //int is = cevent->ids[closest[ip]];
503 Add(z,x,q,q);
504 cog[ip] += x*q;
505 qpl[ip] += q;
506 // cout << ip << " -o- "<<is<<endl;
507 }
508 }
509 // -----------------
510 // add +/- one strip
511 // -----------------
512 for(int ih=0;ih<cevent->GetN(); ih++){
513 float x = cevent->xycoord[ih];
514 float z = cevent->zcoord[ih];
515 float q = cevent->q[ih];
516 int ip = cevent->idp[ih];
517 int is = cevent->ids[ih];
518 if( closest[ip]!=-1 ){
519 int isc = cevent->ids[closest[ip]];
520 if( is == isc+1 || is == isc-1 ){
521 Add(z,x,q,q);
522 cog[ip] += x*q;
523 qpl[ip] += q;
524 // cout << ip << " -+- "<<is<<endl;
525 }
526 }
527 }
528 Fit();
529
530 if(GetN()<3)return 0;
531 for(int ip=0; ip<22; ip++){
532 if(qpl[ip]!=0){
533 cog[ip]=cog[ip]/qpl[ip];
534 // Add(z,cog[ip],cog[ip],0.7);
535 }
536 }
537
538 // cout << "(3)"<<endl;
539 // cout << " n. hits :"<<GetN()<<endl;
540 // cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl;
541 // cout << " Chi2 :"<<GetChi2()<<endl;
542 // cout << " Q axis :"<<GetQaxis()<<endl;
543 // cout << " Q strip :"<<GetQstrip()<<endl;
544 // cout << " ---------------------------"<<endl;
545 return 1;
546 };
547
548
549 int FitAxis( CaloLevel1* calo , Int_t view , Float_t toll ){ return FitAxis(calo,view,toll,0); };
550
551
552
553 //-----------------------------------------------------------------
554 // fit the axis (optimized for interacting patterns?? ...forse)
555 //-----------------------------------------------------------------
556 int FitShower( CaloLevel1* calo , Int_t view , Float_t toll , Bool_t usemechanicalalignement ){
557
558 // cout << "CaloAxis::FitAxis(...)"<<endl;
559
560 Set(calo,view,usemechanicalalignement);
561
562 // ------------------------------
563 // first :
564 // ------------------------------
565 for(int ih=0;ih<cevent->GetN(); ih++){
566 float x = cevent->xycoord[ih];
567 //float z = cevent->zcoord[ih];
568 float q = cevent->q[ih];
569 int ip = cevent->idp[ih];
570 cog[ip] += x*q;
571 qpl[ip] += q;
572 // Add(z,x,q);
573 }
574 for(int ip=0; ip<22; ip++){
575 if(qpl[ip]!=0)cog[ip]=cog[ip]/qpl[ip];
576 // cout << ip << " > "<<qpl[ip]<<" "<<cog[ip]<<endl;
577 }
578 // ----------------
579 // plane of maximum
580 // ----------------
581 float qplmax=1;
582 int ipmax = -1;
583 float dmin[22];
584 int closest[22];
585 for(int ip=0; ip<22; ip++)if(qpl[ip]>qplmax){
586 ipmax=ip;
587 qplmax=qpl[ip];
588 dmin[ip]=1000.;//init
589 closest[ip]=-1;//init
590 }
591 for(int ih=0;ih<cevent->GetN(); ih++){
592 float x = cevent->xycoord[ih];
593 float z = cevent->zcoord[ih];
594 float q = cevent->q[ih];
595 int ip = cevent->idp[ih];
596 // if( (ipmax>3 && ip<=ipmax) || ipmax<=3 )Add(z,x,q);
597 if( ip<=ipmax || ip<=3 ){
598 Add(z,x,q);
599 }
600 }
601 if(GetN()<3)return 0;
602 Fit();
603 // cout << " n. hits :"<<GetN()<<endl;
604 // cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl;
605 // cout << " Chi2 :"<<GetChi2()<<endl;
606 // cout << " Q axis :"<<GetQaxis()<<endl;
607 // cout << " Q strip :"<<GetQstrip()<<endl;
608
609 float P0;
610 float P1;
611 int dof;
612
613 // ---------------------------
614 // second : iteration
615 // ---------------------------
616 int niter=0; // number of iterations
617 int pok[22];
618 // prova
619 // pesi lineari
620 float W22=1.;
621 float W1=10.;
622 float b=(W22-W1)/(22.-1.);// (w22 - w1)/(22-1)
623 float a=W1-b;
624
625 for(int i=0; i<22; i++){
626 pok[i]=0;
627 }
628 do{
629 niter++;
630 float ttoll = toll; //tolerance (cm)
631 if(niter==1)ttoll = 10.*toll;
632 dof = GetN(); //previous fit
633 P0 = par[0];
634 P1 = par[1];
635 Reset(); //fit reset
636 for(int i=0; i<22; i++){
637 cog[i]=0;
638 qpl[i]=0;
639 }
640 for(int ih=0;ih<cevent->GetN(); ih++){
641 float x = cevent->xycoord[ih];
642 float z = cevent->zcoord[ih];
643 float q = cevent->q[ih];
644 int ip = cevent->idp[ih];
645 float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1);
646 // cout<<d<<endl;
647 if(
648 (niter==1 || (niter>1 && pok[ip]==1)) &&
649 (ip<=ipmax || ip<=3) &&
650 // ((ipmax>3 && ip<=ipmax) || ipmax<=3) &&
651 true){
652
653 if(niter==1 && d<dmin[ip]){
654 dmin[ip]=d;
655 closest[ip]=ih;
656 }
657
658 float W=a+b*(ip+1);
659 // cout << ip << " "<<W<<endl;
660 if( d < ttoll || ( niter==2 && ih==closest[ip] )){
661 /* if(ip==ipmax && ipmax>1)Add(z,x,q,W*q); */
662 /* else Add(z,x,q,W); */
663 cog[ip] += x*q;
664 qpl[ip] += q;
665 Add(z,x,q,W);
666 pok[ip]=1;
667 cout << ip << " "<<qpl[ip]<<endl;
668 }
669 }
670 }
671 // break the track if more than 3 planes are missing
672 int nmissing = 0;
673 for(int ip=0; ip<22; ip++){
674 if(pok[ip]==0){
675 nmissing++;
676 }else{
677 nmissing = 0;
678 }
679 if(nmissing==6){
680 for(int ipp=ip+1; ipp<22; ipp++)pok[ipp]=0;
681 break;
682 }
683 }
684
685 for(int ip=0; ip<22; ip++)cout << qpl[ip] << " ";
686 cout << endl;
687
688
689 for(int ip=0; ip<22; ip++)if(qpl[ip]!=0)cog[ip]=cog[ip]/qpl[ip];
690
691
692 if(niter==100)break;
693 if(GetN()<3)break;
694
695 Fit();
696
697 }while(niter==1 || GetN() != dof);
698 // cout << " >>> "<<GetN()<<endl;
699 if(GetN()<3)return 0;
700 Fit();
701 /* // cout << " n. hits :"<<GetN()<<endl; */
702 /* // cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl; */
703 /* // cout << " Chi2 :"<<GetChi2()<<endl; */
704 /* // cout << " Q axis :"<<GetQaxis()<<endl; */
705 /* // cout << " Q strip :"<<GetQstrip()<<endl; */
706
707 /* // --------------------------------------------- */
708 /* // third : selecting only closest strip-clusters */
709 /* // and fit baricenters */
710 /* // --------------------------------------------- */
711 /* P0 = par[0]; */
712 /* P1 = par[1]; */
713 /* Reset(); */
714
715 /* float dmin[22]; */
716 /* int closest[22]; */
717 /* for(int ip=0; ip<22; ip++){ */
718 /* dmin[ip]=999; */
719 /* closest[ip]=-1; */
720 /* } */
721 /* for(int ih=0;ih<cevent->GetN(); ih++){ */
722 /* float x = cevent->xycoord[ih]; */
723 /* float z = cevent->zcoord[ih]; */
724 /* float q = cevent->q[ih]; */
725 /* int ip = cevent->idp[ih]; */
726 /* int is = cevent->ids[ih]; */
727 /* float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1); */
728 /* if( d < toll && d<dmin[ip] && pok[ip]==1 ){ */
729 /* // Add(z,x,q,q); */
730 /* closest[ip]=ih; */
731 /* // cog[ip] += x*q; */
732 /* // qpl[ip] += q; */
733 /* } */
734 /* } */
735 /* for(Int_t ip=0; ip<22; ip++){ */
736 /* if(closest[ip]!=-1){ */
737 /* float x = cevent->xycoord[closest[ip]]; */
738 /* float z = cevent->zcoord[closest[ip]]; */
739 /* float q = cevent->q[closest[ip]]; */
740 /* int is = cevent->ids[closest[ip]]; */
741 /* Add(z,x,q,q); */
742 /* cog[ip] += x*q; */
743 /* qpl[ip] += q; */
744 /* // cout << ip << " -o- "<<is<<endl; */
745 /* } */
746 /* } */
747 /* // add +/- one strip */
748 /* for(int ih=0;ih<cevent->GetN(); ih++){ */
749 /* float x = cevent->xycoord[ih]; */
750 /* float z = cevent->zcoord[ih]; */
751 /* float q = cevent->q[ih]; */
752 /* int ip = cevent->idp[ih]; */
753 /* int is = cevent->ids[ih]; */
754 /* if( closest[ip]!=-1 ){ */
755 /* int isc = cevent->ids[closest[ip]]; */
756 /* if( is == isc+1 || is == isc-1 ){ */
757 /* Add(z,x,q,q); */
758 /* cog[ip] += x*q; */
759 /* qpl[ip] += q; */
760 /* // cout << ip << " -+- "<<is<<endl; */
761 /* } */
762 /* } */
763 /* } */
764 /* Fit(); */
765
766 /* if(GetN()<3)return 0; */
767 /* for(int ip=0; ip<22; ip++){ */
768 /* if(qpl[ip]!=0){ */
769 /* cog[ip]=cog[ip]/qpl[ip]; */
770 /* // Add(z,cog[ip],cog[ip],0.7); */
771 /* } */
772 /* } */
773
774 // cout << " n. hits :"<<GetN()<<endl;
775 // cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl;
776 // cout << " Chi2 :"<<GetChi2()<<endl;
777 // cout << " Q axis :"<<GetQaxis()<<endl;
778 // cout << " Q strip :"<<GetQstrip()<<endl;
779 // cout << " ---------------------------"<<endl;
780 return 1;
781 };
782
783
784 int FitShower( CaloLevel1* calo , Int_t view , Float_t toll ){ return FitAxis(calo,view,toll,0); };
785
786
787 };
788
789 //===============================================================================
790 //
791 //
792 //
793 //
794 //===============================================================================
795
796 #endif // CALOAXIS2_H_

  ViewVC Help
Powered by ViewVC 1.1.23