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

Contents of /calo/flight/CaloFranzini/src/Calib.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.8 - (show annotations) (download)
Fri Jan 25 15:09:07 2008 UTC (16 years, 10 months ago) by mocchiut
Branch: MAIN
Changes since 1.7: +97 -34 lines
it works

1 #include <stdlib.h>
2 #include <iostream>
3 #include <iomanip>
4 //
5 #include <TString.h>
6 #include <TH1F.h>
7 #include <TH2F.h>
8 #include <TMatrixD.h>
9 #include <TMatrixF.h>
10 #include <TArrayF.h>
11 #include <TArrayI.h>
12 //
13 #include <PamLevel2.h>
14 #include <CaloFranzini.h>
15 //
16 using namespace std;
17 //
18 extern Bool_t MATRIX;
19 extern Bool_t FULL;
20 extern Bool_t SIMU;
21 extern Bool_t CRIG;
22 extern Bool_t SRIG;
23 CaloFranzini *cf;
24 Int_t nbin;
25 Float_t rig[18];
26 Float_t rmean[17];
27 Int_t ntot[17];
28 Int_t MDIM = 8213;
29 //Int_t MDIM = 4128;
30 //Float_t qqplane[17][43];
31 //Int_t nnqplane[17][43];
32
33 TArrayF *qplane[17];
34 TArrayI *nqplane[17];
35 TMatrixD *matrix[17];
36 TMatrixD *nmat[17];
37
38 //TMatrixD *fqplane;
39 //TMatrixD *fnqplane;
40 TMatrixD *fqplane[17];
41 TMatrixD *fnqplane[17];
42 TMatrixD *fmatrix[17];
43 //TMatrixF *fnmat[17];
44 //TMatrixD *fmatrix;
45 //TMatrixD *fnmat;
46 //TMatrixF *fmatrix;
47 //TMatrixF *fnmat;
48 TMatrixF *fnmat[17];
49 //Int_t finmat[43][191];
50
51 //===============================================================================
52 bool Select( PamLevel2* event ){
53
54 //---------------------------------------------------------
55 // single track
56 //---------------------------------------------------------
57 if( event->GetTrkLevel2()->GetNTracks()!=1 ) return false;
58 PamTrack *track = event->GetTrack(0);
59 if(!track)return false;
60
61 if ( SIMU ) return true;
62 //------------------------------------------------------------------
63 // tracker pre-selection
64 //------------------------------------------------------------------
65 TrkTrack *trk = track->GetTrkTrack();
66 float rigidity = trk->GetRigidity();
67 if ( CRIG ) rigidity = event->GetCaloLevel2()->qtot/260.;
68 if ( SRIG ) rigidity = event->GetGPamela()->P0;
69 bool TRACK__OK = false;
70 if(
71 trk->chi2 >0 &&
72 trk->GetNX()>=4 &&
73 trk->GetNY()>=3 &&
74 trk->GetLeverArmX()>=5 &&
75 true ) TRACK__OK = true;
76
77 if( !TRACK__OK )return false;
78
79 //------------------------------------------------------------------
80 // TOF pre-selection
81 //------------------------------------------------------------------
82 bool TOF__OK = false;
83 if(
84 event->GetToFLevel2()->GetNHitPaddles(0) == 1 &&
85 event->GetToFLevel2()->GetNHitPaddles(1) == 1 &&
86 event->GetToFLevel2()->GetNHitPaddles(2) == 1 &&
87 event->GetToFLevel2()->GetNHitPaddles(3) == 1 &&
88 event->GetToFLevel2()->GetNHitPaddles(4) >= 1 &&
89 event->GetToFLevel2()->GetNHitPaddles(5) >= 1 &&
90 event->GetToFLevel2()->npmt() <= 18 &&
91 !event->GetAcLevel2()->CARDhit() &&
92 !event->GetAcLevel2()->CAThit() &&
93 true ) TOF__OK = true;
94 if( !TOF__OK && !SIMU)return false;
95 //------------------------------------------------------
96 // no albedo
97 //------------------------------------------------------
98 if( !SIMU && (track->GetToFTrack()->beta[12]<=0.2 ||
99 track->GetToFTrack()->beta[12] >= 1.5) ) return false;
100
101 //------------------------------------------------------
102 bool CUT1 = false;
103 if(
104 trk->nstep<100 &&
105 rigidity<400. &&
106 rigidity>0.1 &&
107 trk->resx[0]<0.001 &&
108 trk->resx[5]<0.001 &&
109 track->IsSolved() &&
110 trk->IsInsideCavity() &&
111 true ) CUT1 = true;
112 //------------------------------------------------------
113 if( !CUT1 )return false;
114 if ( trk->GetDeflection()>0. && !SIMU ) return false;
115
116 //
117 // ELENA'S CUT
118 //
119 //
120 // lever-arm 6
121 //====================================================
122 bool LX6=false;
123 if(
124 track->GetTrkTrack()->GetLeverArmX()==6 &&
125 !track->GetTrkTrack()->IsBad(0,0) &&
126 !track->GetTrkTrack()->IsBad(5,0) &&
127 track->GetTrkTrack()->resx[0]<0.001 &&
128 track->GetTrkTrack()->resx[5]<0.001 &&
129 track->GetTrkTrack()->IsInsideCavity() &&
130 true ) LX6 = true;
131
132 //====================================================
133 // lever-arm 5
134 //====================================================
135 bool LX5=false;
136 if(
137 track->GetTrkTrack()->GetLeverArmX()==5 &&
138 true ){
139 if(
140 track->GetTrkTrack()->XGood(0) && track->GetTrkTrack()->XGood(4)
141 ){
142
143 if(
144 !track->GetTrkTrack()->IsBad(0,0) &&
145 !track->GetTrkTrack()->IsBad(4,0) &&
146 track->GetTrkTrack()->resx[0]<0.001 &&
147 track->GetTrkTrack()->resx[4]<0.001 &&
148 track->GetTrkTrack()->IsInsideCavity() &&
149 true) LX5 = true;
150 }else if (
151 track->GetTrkTrack()->XGood(1) && track->GetTrkTrack()->XGood(5)
152 ){
153
154 if(
155 !track->GetTrkTrack()->IsBad(1,0) &&
156 !track->GetTrkTrack()->IsBad(5,0) &&
157 track->GetTrkTrack()->resx[1]<0.001 &&
158 track->GetTrkTrack()->resx[5]<0.001 &&
159 track->GetTrkTrack()->IsInsideCavity() &&
160 true) LX5 = true;
161 }
162 }
163 if ( !LX5 && !LX6 ) return false;
164 Float_t defl = trk->GetDeflection();
165 float p0 = 1.111588e+00;
166 float p1 = 1.707656e+00;
167 float p2 = 1.489693e-01;
168 float chi2m025 = p0 + fabs(defl)*p1 + defl*defl*p2;
169
170 float def_0 = 0.07;
171 float chi2m025_0 = p0 + fabs(def_0)*p1 + def_0*def_0*p2;
172
173 // int nchi2cut=5;
174 float chi2cut=3.;
175 float chi2m = pow( chi2m025-chi2m025_0+pow(chi2cut,0.25), 4.);
176 bool CUT2 = false;
177 if(
178 track->GetTrkTrack()->chi2 < chi2m &&
179 true ) CUT2 = true;
180 if ( !CUT2 ) return false;
181 float dedxtrk = trk->GetDEDX();
182 // float zcutn = 9. + 20./(rigidity*rigidity);
183 float zcut2 = 3. + 4.3/(rigidity*rigidity);
184 float zcut1 = 0.52 + 0.455/(rigidity*rigidity);
185 Bool_t Z1 = false;
186 if(dedxtrk > zcut1 && dedxtrk < zcut2){
187 Z1=true;
188 }
189 if ( !Z1 && !SIMU ) return false;
190 //------------------------------------------------------
191 //
192 // energy momentum match
193 //
194 Float_t qtotimp = event->GetCaloLevel2()->qtot / trk->GetRigidity();
195 Float_t qcut2 = (-0.5 * trk->GetRigidity() + 150.) * 1.1;
196 if ( qcut2 < 55. ) qcut2 = 55.;
197 if ( qtotimp <= qcut2 ) return false;
198 //
199 for (Int_t i=0; i < 22; i++){
200 if ( track->GetCaloTrack()->tibar[i][1] < 0 || track->GetCaloTrack()->tibar[i][0] < 0 ){
201 return false;
202 };
203 };
204 //
205 if ( event->GetCaloLevel2()->qtot == 0. ) return false;
206 if ( rigidity>5. && track->GetCaloTrack()->qtrack/event->GetCaloLevel2()->qtot < 0.4 ) return false;
207 if ( rigidity<1. && track->GetToFTrack()->beta[12] < 0.8 ) return false;
208 if ( rigidity>50. ){
209 if ( trk->GetNX()<5 &&
210 trk->GetNY()<4 ) return false;
211 //
212 Bool_t sphit = false;
213 for ( Int_t plane = 0; plane < 6; plane++){
214 if ( !trk->XGood(plane) ){
215 for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsx(); sing++){
216 TClonesArray &t = *(event->GetTrkLevel2()->SingletX);
217 TrkSinglet *singlet = (TrkSinglet*)t[sing];
218 if ( (singlet->plane-1) == plane ){
219 Float_t x = (singlet->coord[0]+singlet->coord[1])/2.;
220 if ( fabs(track->GetTrkTrack()->xv[plane] - x) < 1. ) sphit = true;
221 };
222 };
223 };
224 if ( !trk->YGood(plane) ){
225 for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsy(); sing++){
226 TClonesArray &t = *(event->GetTrkLevel2()->SingletY);
227 TrkSinglet *singlet = (TrkSinglet*)t[sing];
228 if ( (singlet->plane-1) == plane ){
229 Float_t x1 = (singlet->coord[0]);
230 Float_t x2 = (singlet->coord[1]);
231 if ( fabs(track->GetTrkTrack()->yv[plane] - x1) < 1. ) sphit = true;
232 if ( fabs(track->GetTrkTrack()->yv[plane] - x2) < 1. ) sphit = true;
233 };
234 };
235 };
236 };
237 if ( sphit ) return false; // spurious hit along the track
238 };
239 //
240 Int_t ti0 = track->GetCaloTrack()->tibar[0][1]-1;
241
242 Int_t view = 0;
243 Int_t plane = 0;
244 Int_t strip = 0;
245 Float_t mip = 0.;
246 //
247 for ( Int_t i=0; i<event->GetCaloLevel1()->istrip; i++ ){
248 //
249 mip = event->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
250 if ( view == 1 && plane == 0 && strip == ti0 && mip > 4.) return false;
251 if ( view == 1 && (plane >0 || strip > ti0) ) break;
252 };
253 // if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;
254
255 // if ( rigidity > 2.2 || rigidity < 1.5 ) return false;
256 // printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG);
257 //
258 return true;
259 }
260 //===============================================================
261 // Create histograms
262 //
263 //
264 //
265 //
266 //
267 //===============================================================
268 void CreateHistos( PamLevel2* event , TString file){
269
270 cf = new CaloFranzini(event);
271 //
272 if ( FULL ){
273 cf->CalculateLongTZeta(false);
274 cf->CalculateFullTZeta();
275 };
276 //
277 if ( MATRIX ){
278 cf->UpdateMatrixFile(file.Data());
279 if ( !FULL ){
280 cf->LoadBin();
281 cf->LoadLong();
282 } else {
283 cf->LoadBin(true);
284 cf->LoadFull();
285 };
286 } else {
287 cf->CreateMatrixFile(file.Data());
288 };
289 //
290 //
291 nbin = 18;
292 rig[0] = 0.1;
293 rig[1] = 0.5;
294 rig[2] = 1.;
295 rig[3] = 1.5;
296 rig[4] = 2.2;
297 rig[5] = 3.;
298 rig[6] = 4.;
299 rig[7] = 5.;
300 rig[8] = 6.;
301 rig[9] = 8.;
302 rig[10] = 10.;
303 rig[11] = 15.;
304 rig[12] = 25.;
305 rig[13] = 35.;
306 rig[14] = 50.;
307 rig[15] = 100.;
308 rig[16] = 200.;
309 rig[17] = 400.;
310 //
311 memset(rmean, 0, 17*sizeof(Float_t));
312 memset(ntot, 0, 17*sizeof(Int_t));
313 // memset(finmat, 0, 43*191*sizeof(Int_t));
314 // Double_t tol = 1E-20;
315 //
316 for (Int_t i=0; i < 17 ; i++){
317 // for (Int_t i=3; i < 4 ; i++){
318 if ( !FULL ){
319 matrix[i] = new TMatrixD(43,43);
320 qplane[i] = new TArrayF(43);
321 nqplane[i] = new TArrayI(43);
322 nmat[i] = new TMatrixD(43,43);
323 } else {
324 if ( MATRIX ){
325 // fmatrix = new TMatrixF(4128,4128);
326 // fnmat = new TMatrixF(4128,4128);
327 // fmatrix = new TMatrixF(8213,8213);
328 // fnmat = new TMatrixF(8213,8213);
329 // fmatrix = new TMatrixF(MDIM,MDIM);
330 // fnmat = new TMatrixF(MDIM,MDIM);
331 // fmatrix[i] = new TMatrixF(1849,1849);
332 // fnmat[i] = new TMatrixF(43,43);
333 // fmatrix[i] = new TMatrixD(1333,1333);
334 // fnmat[i] = new TMatrixF(43,31);
335 fmatrix[i] = new TMatrixD(473,473);
336 fnmat[i] = new TMatrixF(43,11);
337 // fmatrix[i]->SetTol(tol);
338 // cf->WriteFullMatrix(fmatrix, i);
339 // cf->WriteFullNMatrix(fnmat, i);
340 // delete fmatrix;
341 // delete fnmat;
342 //fnmat[i] = new TMatrixI(8213,8213);
343 } else {
344 // fqplane = new TMatrixD(43,191); // 43 planes x 191 strip (= 1 + 95 x 2, one strip is the one transversed by the track that could be on the extreme right or left)
345 // fnqplane = new TMatrixD(43,191);//
346 // fqplane[i] = new TMatrixD(43,43); // 43 planes x 43 "strip", where 43 = 50 + 14 + 5 + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + [1] + ...
347 // fnqplane[i] = new TMatrixD(43,43);// 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
348 //
349 // fqplane[i] = new TMatrixD(43,31); // 43 planes x 43 "strip", where 43 = 50 + 14 + 6 + 5 + 3 + 3 + 3 + 2 + 2 + 2 + 1 + 1 + 1 + 1 + 1 + [1] + 1 + 1 + 1 + 1 + ...
350 // fnqplane[i] = new TMatrixD(43,31);// 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ...
351 fqplane[i] = new TMatrixD(43,11); // 43 planes x 11 "strip", where 43 = 84 + 6 + 3 + 1 + 1 +[1]+ 1 + 1 + 3 + 6 + 84
352 fnqplane[i] = new TMatrixD(43,11);// 0 1 2 3 4 5 6 7 8 9 10
353 //
354 // cf->WriteFullMean(fqplane, i);
355 // cf->WriteFullNMean(fnqplane, i);
356 // delete fqplane;
357 // delete fnqplane;
358 //
359 };
360 };
361 };
362 //
363 }
364
365 //===============================================================
366 void FindAverage( PamLevel2* L2, int iev ){
367 //
368 Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
369 if ( SRIG ) erig = L2->GetGPamela()->P0;
370 if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.;
371 //
372 Int_t rbi = 0;
373 for (Int_t i = 0; i<nbin-1; i++){
374 if ( erig>=rig[i] && erig < rig[i+1] ){
375 rbi = i;
376 break;
377 };
378 };
379 //
380 if ( erig < rig[0] ) return;
381 if ( erig >= rig[nbin-1] ) return;
382 //
383 rmean[rbi] += erig;
384 ntot[rbi]++;
385 //
386 if (!FULL ){
387 Int_t dgf = 43;
388 //
389 for (Int_t i=0; i<dgf; i++){
390 (*nqplane[rbi])[i]++;
391 };
392 //
393 // Fill the estrip matrix
394 //
395 Int_t nplane = 0;
396 Int_t view = 0;
397 Int_t plane = 0;
398 Int_t strip = 0;
399 Float_t mip = 0.;
400 //
401 for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
402 //
403 mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
404 //
405 nplane = 1 - view + 2 * plane;
406 if ( erig > 4. && nplane == 0 && mip > 15. ) printf(" IEV %i erig %f OBT %u pkt %u file %s \n",iev,erig,L2->GetOrbitalInfo()->OBT,L2->GetOrbitalInfo()->pkt_num,L2->GetPamTree()->GetFile()->GetName());
407 //printf(" IEV %i OBT %u pkt %u file %s \n",iev,L2->GetOrbitalInfo()->OBT,L2->GetOrbitalInfo()->pkt_num,L2->GetPamTree()->GetFile()->GetName());
408 if ( nplane > 37 ) nplane--;
409 if ( nplane < dgf ){
410 (*qplane[rbi])[nplane] += mip;
411 };
412 //
413 };
414 } else {
415 //
416 // FULL CALORIMETER
417 //
418 // fqplane = cf->LoadFullAverage(rbi);
419 // fnqplane = cf->LoadFullNAverage(rbi);
420 CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
421 //
422 Int_t nplane = 0;
423 Int_t view = 0;
424 Int_t plane = 0;
425 Int_t strip = 0;
426 Float_t mip = 0.;
427 //
428 Int_t cs = 0;
429 Int_t cd = 0;
430 Int_t mstrip = 0;
431 //
432 for (Int_t j=0; j<2; j++){
433 for (Int_t i=0; i<22; i++){
434 nplane = 1 - j + 2*i;
435 if ( nplane > 37 ) nplane--;
436 //
437 cs = ct->tibar[i][j] - 1;
438 //
439 cd = 95 - cs;
440 //
441 Int_t oldstr = -1;
442 for (Int_t k=0; k<191; k++){
443 mstrip = cd + k;
444 // if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.;
445 // if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.;
446 Int_t lstr = cf->ConvertStrip(mstrip);
447 if ( oldstr != lstr ){
448 (*fnqplane[rbi])[nplane][lstr] += 1.;
449 oldstr = lstr;
450 };
451 };
452 };
453 };
454 //
455 //
456 for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
457 //
458 mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
459 //
460 nplane = 1 - view + 2 * plane;
461 if ( nplane > 37 ) nplane--;
462 //
463 cs = ct->tibar[plane][view] - 1;
464 //
465 cd = 95 - cs;
466 //
467 mstrip = cd + strip;
468 //
469 Int_t lstr = cf->ConvertStrip(mstrip);
470 // (*fqplane[rbi])[nplane][mstrip] += mip;
471 // (*fqplane)[nplane][mstrip] += mip;
472 (*fqplane[rbi])[nplane][lstr] += mip;
473 //
474 };
475 //
476 // cf->WriteFullMean(fqplane, rbi);
477 // cf->WriteFullNMean(fnqplane, rbi);
478 // cf->UnLoadFullAverage(rbi);
479 // cf->UnLoadFullNAverage(rbi);
480 // delete fqplane;
481 // delete fnqplane;
482 //
483 };
484 }
485
486 void CalculateAverage(){
487 //
488 if ( !FULL ){
489 for (Int_t i=0; i<nbin-1; i++){
490 if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];
491 for (Int_t j=0; j<43 ; j++){
492 if ( (*nqplane[i])[j] > 0 ){
493 (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];
494 } else {
495 (*qplane[i])[j] = 0.;
496 };
497 printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]);
498 };
499 };
500 for (Int_t i=0; i<nbin-1; i++){
501 //
502 cf->WriteLongMean(qplane[i], i);
503 //
504 };
505 } else {
506 //
507 for (Int_t i=0; i<nbin-1; i++){
508 // fqplane = cf->LoadFullAverage(i);
509 // fnqplane = cf->LoadFullNAverage(i);
510 if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]);
511 //
512 for (Int_t j=0; j<43 ; j++){
513 // for (Int_t k=0; k<191; k++){
514 // for (Int_t k=0; k<43; k++){
515 // for (Int_t k=0; k<31; k++){
516 for (Int_t k=0; k<11; k++){
517 // if ( (*fnqplane[i])[j][k] > 0 ){
518 // (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
519 // } else {
520 // (*fqplane[i])[j][k] = 0.;
521 // };
522 // printf(" BIN %i plane %i strip %i average energy %f qplane %f nqplane %f \n",i,j,k,rmean[i],(*fqplane[i])[j][k],(*fnqplane[i])[j][k]);
523 if ( (*fnqplane[i])[j][k] > 0 ){
524 if ( (*fqplane[i])[j][k] == 0. ) (*fqplane[i])[j][k] = 0.7;
525 (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
526 } else {
527 (*fqplane[i])[j][k] = 0.;
528 };
529 // printf(" BIN %i plane %i strip %i average energy %f qplane %f nqplane %f \n",i,j,k,rmean[i],(*fqplane)[j][k],(*fnqplane)[j][k]);
530 };
531 };
532 cf->WriteFullMean(fqplane[i], i);
533 cf->WriteFullNMean(fnqplane[i], i);
534 // cf->UnLoadFullAverage(i);
535 // cf->UnLoadFullNAverage(i);
536 // delete fqplane;
537 // delete fnqplane;
538 };
539 //
540 // for (Int_t i=0; i<nbin-1; i++){
541 // //
542 // cf->WriteFullMean(fqplane[i], i);
543 // //
544 // };
545 };
546 //
547 cf->WriteNumBin(nbin);
548 //
549 TArrayF *rigbin = new TArrayF(18, rig);
550 cf->WriteRigBin(rigbin);
551 //
552 TArrayF *rmeanbin = new TArrayF(17, rmean);
553 if ( FULL ){
554 TFile *file = cf->GetFullFile();
555 file->cd();
556 file->WriteObject(&(*rmeanbin), "binrigmean");
557 } else {
558 TFile *file = cf->GetFile();
559 file->cd();
560 file->WriteObject(&(*rmeanbin), "binrigmean");
561 };
562 //
563 //
564 }
565
566 //===============================================================
567 void FindMatrix( PamLevel2* L2, int iev ){
568 //
569 Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
570 if ( SRIG ) erig = L2->GetGPamela()->P0;
571 if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.;
572 //
573 Int_t rbi = 0;
574 for (Int_t i = 0; i<nbin-1; i++){
575 if ( erig>=rig[i] && erig < rig[i+1] ){
576 rbi = i;
577 break;
578 };
579 };
580 //
581 if ( erig < rig[0] ) return;
582 if ( erig >= rig[nbin-1] ) return;
583 //
584 if ( !FULL ){
585 Int_t dgf = 43;
586 //
587 for (Int_t i=0; i<dgf; i++){
588 for (Int_t j=0; j<dgf; j++){
589 (*nmat[rbi])[i][j] += 1.;
590 };
591 };
592 //
593 // Fill the estrip matrix
594 //
595 Int_t nplane = 0;
596 Int_t view = 0;
597 Int_t plane = 0;
598 Int_t strip = 0;
599 Float_t mip = 0.;
600 Float_t hpl[43];
601 memset(hpl,0,43*sizeof(Float_t));
602 for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
603 //
604 mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
605 //
606 nplane = 1 - view + 2 * plane;
607 if ( nplane > 37 ) nplane--;
608 if ( nplane < dgf ){
609 hpl[nplane] += mip;
610 };
611 //
612 };
613 //
614 for (Int_t i=0; i<dgf; i++){
615 for (Int_t j=0; j<dgf; j++){
616 (*matrix[rbi])[i][j] += (hpl[i] - cf->GetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig));
617 };
618 };
619 } else {
620 //
621 // FULL CALORIMETER
622 //
623 // if ( rbi != 3 ) return;
624 printf(" matrix %i IEV %i \n",rbi,iev);
625 // fmatrix = cf->LoadFullMatrix(rbi);
626 // cf->LoadFullMatrix(rbi,fmatrix);
627 // fnmat = cf->LoadFullNMatrix(rbi);
628 // printf(" done \n");
629 // printf(" start loop \n");
630 //
631 CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
632 //
633 Int_t nplane = 0;
634 Int_t view = 0;
635 Int_t plane = 0;
636 Int_t strip = 0;
637 Float_t mip = 0.;
638 //
639 // Int_t mindgf = 48;
640 // Int_t dgf = 143;
641 // Int_t mindgf = 0; //tutto
642 // Int_t dgf = 191; //tutto
643 // Int_t mindgf = 94;
644 // Int_t dgf = 96;
645 // Int_t mindgf = 84;
646 // Int_t dgf = 106;
647 Int_t mindgf = 0;
648 // Int_t dgf = 43;
649 Int_t dgf = 11;
650 Int_t cs = 0;
651 Int_t cd = 0;
652 Int_t mstrip = 0;
653 //
654 // Float_t mipv[43][43];
655 // memset(mipv,0,43*43*sizeof(Float_t));
656 // Float_t mipv[43][31];
657 // memset(mipv,0,43*31*sizeof(Float_t));
658 Float_t mipv[43][11];
659 memset(mipv,0,43*11*sizeof(Float_t));
660 //
661 for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
662 //
663 mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
664 //
665 nplane = 1 - view + 2 * plane;
666 if ( nplane > 37 ) nplane--;
667 //
668 cs = ct->tibar[plane][view] - 1;
669 //
670 cd = 95 - cs;
671 //
672 mstrip = cd + strip;
673 //
674 Int_t lstr = cf->ConvertStrip(mstrip);
675 mipv[nplane][lstr] += mip;
676 //
677 };
678 //
679 Float_t mip1 = 1.;
680 Int_t cs1;
681 Int_t cd1;
682 Float_t mip2 = 1.;
683 Int_t cs2;
684 Int_t cd2;
685 Int_t mi = -1;
686 Int_t mj = -1;
687 Int_t nn1 = 0;
688 Int_t pl1 = 0;
689 Int_t vi1 = 0;
690 Int_t nn2 = 0;
691 Int_t pl2 = 0;
692 Int_t vi2 = 0;
693 Int_t mstrip1min = 0;
694 Int_t mstrip1max = 0;
695 Int_t mstrip2min = 0;
696 Int_t mstrip2max = 0;
697 //
698 Int_t toto = 0;
699 //
700 //
701 //
702 Int_t therigb = rbi;
703 //
704 if ( erig < rig[0] ){
705 therigb = 0;
706 };
707 if ( erig >= rig[nbin-2] ){
708 therigb = nbin - 3;
709 };
710 //
711 for (Int_t nplane1=0; nplane1<43; nplane1++){
712 if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
713 vi1 = 1;
714 if ( nn1%2 ) vi1 = 0;
715 pl1 = (nn1 - 1 + vi1)/2;
716 //
717 cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1
718 //
719 cd1 = 95 - cs1;
720 //
721 Int_t at1 = TMath::Max(0,(cd1+0));
722 Int_t at2 = TMath::Min(190,(cd1+95));
723 mstrip1min = cf->ConvertStrip(at1);
724 mstrip1max = cf->ConvertStrip(at2) + 1;
725 // mstrip1min = cf->ConvertStrip(TMath::Max(mindgf,(cd1+0)));
726 // mstrip1max = cf->ConvertStrip(TMath::Min(dgf,(cd1+95))) + 1;
727 //
728 if ( nplane1 == 0 || nplane1 == 42 ) printf(" pl %i mstrip1min %i mstrip1max %i mindgf %i dgf %i cd1 %i\n",nplane1,mstrip1min,mstrip1max,mindgf,dgf,cd1);
729 //
730 for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
731 // printf(".\n");
732 //
733 mj = -1;
734 //
735 mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,therigb);
736 //
737 // mi = (nplane1 * 191) + mstrip1;
738 // mi = (nplane1 * 43) + mstrip1;
739 // mi = (nplane1 * 31) + mstrip1;
740 mi = (nplane1 * 11) + mstrip1;
741 //
742 // if ( mstrip1 > mstrip1min ) break;
743 // if ( mstrip1 > dgf ) break;
744 // if ( mstrip1 >= mindgf && mstrip1 <= dgf && mstrip1 >= mstrip1min && mstrip1 <= mstrip1max ){
745 //
746 // finmat[nplane1][mstrip1]++;
747 (*fnmat[rbi])[nplane1][mstrip1] += 1.;
748 //
749 if ( mip1 != 0. ){
750 //
751 for (Int_t nplane2=0; nplane2<43; nplane2++){
752 //
753 if ( nplane2 >= 37 ) nn2 = nplane2 + 1;
754 vi2 = 1;
755 if ( nn2%2 ) vi2 = 0;
756 pl1 = (nn2 - 1 + vi2)/2;
757 //
758 cs2 = ct->tibar[pl2][vi2] - 1;
759 //
760 cd2 = 95 - cs2;
761 //
762 // mstrip2min = cd2 + 0;
763 // mstrip2max = cd2 + 95;
764 Int_t t1 = TMath::Max(0,(cd2+0));
765 Int_t t2 = TMath::Min(190,(cd2+95));
766 mstrip2min = cf->ConvertStrip(t1);
767 mstrip2max = cf->ConvertStrip(t2) + 1;
768 //
769 if ( nplane1 == 0 && nplane2 == 0 && mstrip1==mstrip1min ) printf(" mstrip2min %i mstrip2max %i \n",mstrip2min,mstrip2max);
770 //
771 for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
772 //
773 mip2 = mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig,therigb);
774 //
775 if ( mip2 != 0. ){
776 //
777 // mj = (nplane2 * 191) + mstrip2;
778 // mj = (nplane2 * 43) + mstrip2;
779 // mj = (nplane2 * 31) + mstrip2;
780 // Int_t sh = -15 + nplane1;
781 // if ( sh > 15 ) sh -= 31*nplane1;
782 // //
783 // mj = (nplane2 * 31) + mstrip2 + sh;
784 // //
785 // if ( mj < 0 ) mj += 1333;
786 // if ( mj >= 1333 ) mj -= 1333;
787 //
788 //
789 //
790 //
791 // Int_t sh = -5 + nplane1;
792 // if ( sh > 5 ) sh -= 11*nplane1;
793 // //
794 // mj = (nplane2 * 11) + mstrip2 + sh;
795 // //
796 // if ( mj < 0 ) mj += 473;
797 // if ( mj >= 473 ) mj -= 473;
798 // //
799 //
800 mj = (nplane2 * 11) + mstrip2;
801 //
802 // printf(" mi %i mj %i sh %i \n",mi,mj,sh);
803 //
804 // mj++;
805 //
806 // if ( mstrip2 > mstrip2min ) break;
807 // if ( mstrip2 > dgf ) break;
808 // if ( mstrip2 >= mindgf && mstrip2 <= dgf && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max ){
809 // if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){
810 // (*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig));
811 // (*fnmat[rbi])[mi][mj] += 1.;
812 (*fmatrix[rbi])[mi][mj] += (mip1 * mip2); // giusto
813 // (*fmatrix)[mi][mj] += (mip1 * mip2) * 1000000.;
814 toto++;
815 // (*fmatrix)[mi][mj] += 1.;
816 // cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);
817 // cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi);
818 // (*fnmat)[mi][mj] += 1.;
819 // };
820 };
821 };
822 };
823 };
824 };
825 };
826 //
827 printf(" toto = %i \n",toto);
828 printf("\n done \n");
829 // printf(" write matrix \n");
830 // cf->WriteFullMatrix(fmatrix, rbi);
831 // cf->WriteFullNMatrix(fnmat, rbi);
832 // printf(" done \n");
833 // printf(" unload matrix \n");
834 // cf->UnLoadFullMatrix(rbi);
835 // cf->UnLoadFullNMatrix(rbi);
836 // printf(" done \n");
837 // printf(" delete matrix \n");
838 // delete fmatrix;
839 // delete fnmat;
840 // printf(" done \n");
841 };
842 }
843
844 //===============================================================
845 // Save histograms
846 //
847 //
848 //
849 //
850 //
851 //===============================================================
852 void SaveHistos(){
853 //
854 if ( MATRIX ){
855 //
856 printf("Finished, calculating average and inverting matrices\n");
857 //
858 if ( !FULL ){
859 for (Int_t i=0; i<nbin-1; i++){
860 //
861 // determine the average matrix
862 //
863 for (Int_t ii=0; ii<43; ii++){
864 for (Int_t j=0; j<43; j++){
865 if ( (*nmat[i])[ii][j] > 0. ){
866 (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];
867 } else {
868 (*matrix[i])[ii][j] = 0.;
869 };
870 };
871 };
872 //
873 cf->WriteLongMatrix(matrix[i],i);
874 //
875 if ( matrix[i]->Determinant() == 0. ){
876 printf("\n");
877 for (Int_t ii=0; ii<43; ii++){
878 for (Int_t j=0; j<43; j++){
879 printf(" %.f",(*matrix[i])[ii][j]);
880 };
881 printf("\n");
882 };
883 printf("\n");
884 printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
885 } else {
886 Double_t det = 0.;
887 TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));
888 printf(" Bin %i determinant is %f \n",i,det);
889 cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);
890 };
891 };
892 } else {
893 //
894 // FULL
895 //
896 for (Int_t i=0; i<nbin-1; i++){
897 // for (Int_t i=3; i<5; i++){
898 // for (Int_t i=3; i<4; i++){
899 //
900 // determine the average matrix
901 //
902 // fmatrix = cf->LoadFullMatrix(i);
903 // fnmat = cf->LoadFullNMatrix(i);
904 //
905 // for (Int_t ii=0; ii<MDIM; ii++){
906 // for (Int_t j=0; j<MDIM; j++){
907 // // if ( (*fnmat[i])[ii][j] > 0. ){
908 // // (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
909 // // } else {
910 // // (*fmatrix[i])[ii][j] = 0.;
911 // // };
912 // if ( (*fnmat)[ii][j] > 0. ){
913 // (*fmatrix)[ii][j] /= (*fnmat)[ii][j];
914 // } else {
915 // (*fmatrix)[ii][j] = 0.;
916 // };
917 // };
918 // };
919 //
920 // TMatrixD *mymat3 = new TMatrixD(129,129);
921 // TMatrixD *mymat5 = new TMatrixD(215,215);
922 // TMatrixD *mymat7 = new TMatrixD(301,301);
923 // TMatrixD *mymat9 = new TMatrixD(387,387);
924 // TMatrixD *mymat11 = new TMatrixD(473,473);
925 // TMatrixD *mymat17 = new TMatrixD(731,731);
926 // TMatrixF *mymat = new TMatrixF(129,129);
927 // TMatrixF *mymat = new TMatrixF(989,989);
928 Int_t i1 = -1;
929 Int_t j1 = -1;
930 // int mi,mj;
931 Int_t nonzero = 0;
932 Int_t nonzero1 = 0;
933 for (Int_t ii=0; ii<43; ii++){
934 // for (Int_t j=0; j<191; j++){
935 // for (Int_t j=0; j<43; j++){
936 // for (Int_t j=0; j<31; j++){
937 for (Int_t j=0; j<11; j++){
938 // if ( (*fnmat[i])[ii][j] > 0. ){
939 // (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
940 // } else {
941 // (*fmatrix[i])[ii][j] = 0.;
942 // };
943 // i1 = (ii * 191) + j;
944 // i1 = (ii * 43) + j;
945 // i1 = (ii * 31) + j;
946 i1 = (ii * 11) + j;
947 // j1 = -1;
948 for (Int_t iij=0; iij<43; iij++){
949 // for (Int_t jj=0; jj<191; jj++){
950 // for (Int_t jj=0; jj<43; jj++){
951 // for (Int_t jj=0; jj<31; jj++){
952 for (Int_t jj=0; jj<11; jj++){
953 //
954 // j1 = (iij * 191) + jj;
955 // j1 = (iij * 43) + jj;
956 // Int_t sh = -15 + ii;
957 // if ( sh > 15 ) sh -= 31*ii;
958 // //
959 // j1 = (iij * 31) + jj + sh;
960 // //
961 // if ( j1 < 0 ) j1 += 1333;
962 // if ( j1 >= 1333 ) j1 -= 1333;
963 //
964 //
965 //
966 // Int_t sh = -5 + ii;
967 // if ( sh > 5 ) sh -= 11*ii;
968 // //
969 // j1 = (iij * 11) + jj + sh;
970 // //
971 // if ( j1 < 0 ) j1 += 473;
972 // if ( j1 >= 473 ) j1 -= 473;
973 //
974 j1 = (iij * 11) + jj;
975 // j1++;
976 // if ( finmat[ii][j] > 0 ){
977 // (*fmatrix)[i1][j1] /= finmat[ii][j];
978 if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix[i])[i1][j1] == 0. || !((*fmatrix[i])[i1][j1] == (*fmatrix[i])[i1][j1]) ){
979 (*fmatrix[i])[i1][j1] = 1.;
980 } else {
981 (*fmatrix[i])[i1][j1] /= (*fnmat[i])[ii][j];
982 nonzero++;
983 if ( i1 == 0 ) nonzero1++;
984 };
985 //
986 // if ( j>=7 && j <=23 && jj >=7 && jj<=23 ){
987 // Int_t mi17 = (ii*3) + j -7;
988 // Int_t mj17 = (iij*3) + jj -7;
989 // (*mymat17)[mi17][mj17] = (*fmatrix[i])[i1][j1];
990 // };
991 // if ( j>=10 && j <=20 && jj >=10 && jj<=20 ){
992 // Int_t mi11 = (ii*3) + j -10;
993 // Int_t mj11 = (iij*3) + jj -10;
994 // (*mymat11)[mi11][mj11] = (*fmatrix[i])[i1][j1];
995 // };
996 // if ( j>=11 && j <=19 && jj >=11 && jj<=19 ){
997 // Int_t mi9 = (ii*3) + j -11;
998 // Int_t mj9 = (iij*3) + jj -11;
999 // (*mymat9)[mi9][mj9] = (*fmatrix[i])[i1][j1];
1000 // };
1001 // if ( j>=12 && j <=18 && jj >=12 && jj<=18 ){
1002 // Int_t mi7 = (ii*3) + j -12;
1003 // Int_t mj7 = (iij*3) + jj -12;
1004 // (*mymat7)[mi7][mj7] = (*fmatrix[i])[i1][j1];
1005 // };
1006 // if ( j>=13 && j <=17 && jj >=13 && jj<=17 ){
1007 // Int_t mi5 = (ii*3) + j -13;
1008 // Int_t mj5 = (iij*3) + jj -13;
1009 // (*mymat5)[mi5][mj5] = (*fmatrix[i])[i1][j1];
1010 // };
1011 // if ( j>=14 && j <=16 && jj >=14 && jj<=16 ){
1012 // Int_t mi3 = (ii*3) + j -14;
1013 // Int_t mj3 = (iij*3) + jj -14;
1014 // (*mymat3)[mi3][mj3] = (*fmatrix[i])[i1][j1];
1015 // };
1016
1017
1018 // if ( j>=94 && j <=96 && jj >=94 && jj<=96 ){
1019 // mi = (ii*3) + j -94;
1020 // mj = (iij*3) + jj -94;
1021 // (*mymat)[mi][mj] = (*fmatrix)[i1][j1];
1022 // };
1023
1024
1025 // if ( j>=84 && j <=106 && jj >=84 && jj<=106 ){
1026 // mi = (ii*3) + j -84;
1027 // mj = (iij*3) + jj -84;
1028 // (*mymat)[mi][mj] = (*fmatrix)[i1][j1];
1029 // };
1030
1031 };
1032 };
1033 };
1034 };
1035 //
1036 printf(" Matrix has %i non-zero elements \n",nonzero);
1037 // printf(" Matrix has %i non-zero elements on the first row\n",nonzero1);
1038 //
1039 // Bool_t BAD = false;
1040 // for (Int_t ii=0; ii<43; ii++){
1041 // for (Int_t j=0; j<191; j++){
1042 // //
1043 // i1 = (ii * 191) + j;
1044 // //
1045 // for (Int_t iij=0; iij<43; iij++){
1046 // for (Int_t jj=0; jj<191; jj++){
1047 // //
1048 // j1 = (iij * 191) + jj;
1049 // //
1050 // // printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
1051 // if ( (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1]==(*fmatrix)[i1][j1]) ){
1052 // printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
1053 // printf(" che schifo! \n");
1054 // BAD = true;
1055 // };
1056 // //
1057 // };
1058 // };
1059 // };
1060 // };
1061 // //
1062 // if ( BAD ) printf(" questa matrice fa cagare \n");
1063 //
1064 //
1065 // cf->WriteFullMatrix(fmatrix, i);
1066 // cf->WriteFullNMatrix(fnmat, i);
1067 //cf->WriteFullMatrix(fmatrix[i],i);// OK
1068 // cf->WriteFullNMatrix(fnmat[i], i); //OK
1069 //
1070 // TDecompSVD svd(*fmatrix[i]);
1071 // Bool_t ok = svd.Decompose();
1072 //
1073 Double_t zero = (Double_t)0.0;
1074 //
1075 if ( fmatrix[i]->Determinant() == zero ){
1076 //if ( fmatrix->Determinant() == 0. ){
1077 printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1078 cf->WriteFullMatrix(fmatrix[i],i);// OK
1079 cf->WriteFullNMatrix(fnmat[i], i); //OK
1080 } else {
1081 // };
1082 // if ( i == 3 ){
1083 // if ( ok ){
1084 // Double_t tol = 1E-20;
1085 // TDecompSVD svd((*fmatrix)[i],tol);
1086 // svd.Decompose();
1087 // TMatrixD svdInv = svd.Invert();
1088 // svdInv.Print("svdInv");
1089 // cout << "condition: " << svd.Condition() << endl;
1090 // cf->WriteInvertedFullMatrix((TMatrixD)svdInv,999);
1091
1092 Double_t det = 0.;
1093 TMatrixD invmatrix = (TMatrixD)(fmatrix[i]->Invert(&det));
1094 printf(" Bin %i determinant is %f \n",i,det);
1095 cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,i);
1096 };
1097
1098 // if ( mymat3->Determinant() == 0. ){
1099 // printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1100 // } else {
1101 // Double_t det = 0.;
1102 // TMatrixD invmatrix = (TMatrixD)(mymat3->Invert(&det));
1103 // printf(" Mymat3 determinant is %f \n",det);
1104 // cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1103);
1105 // };
1106 // cf->WriteFullMatrix(mymat3, 103);
1107 // if ( mymat5->Determinant() == 0. ){
1108 // printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1109 // } else {
1110 // Double_t det = 0.;
1111 // TMatrixD invmatrix = (TMatrixD)(mymat5->Invert(&det));
1112 // printf(" Mymat5 determinant is %f \n",det);
1113 // cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1105);
1114 // };
1115 // cf->WriteFullMatrix(mymat5, 105);
1116 // if ( mymat7->Determinant() == 0. ){
1117 // printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1118 // } else {
1119 // Double_t det = 0.;
1120 // TMatrixD invmatrix = (TMatrixD)(mymat7->Invert(&det));
1121 // printf(" Mymat7 determinant is %f \n",det);
1122 // cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1107);
1123 // };
1124 // cf->WriteFullMatrix(mymat7, 107);
1125 // if ( mymat9->Determinant() == 0. ){
1126 // printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1127 // } else {
1128 // Double_t det = 0.;
1129 // TMatrixD invmatrix = (TMatrixD)(mymat9->Invert(&det));
1130 // printf(" Mymat3 determinant is %f \n",det);
1131 // cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1109);
1132 // };
1133 // cf->WriteFullMatrix(mymat9, 109);
1134 // if ( mymat11->Determinant() == 0. ){
1135 // printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1136 // } else {
1137 // Double_t det = 0.;
1138 // TMatrixD invmatrix = (TMatrixD)(mymat11->Invert(&det));
1139 // printf(" Mymat11 determinant is %f \n",det);
1140 // cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1111);
1141 // };
1142 // cf->WriteFullMatrix(mymat11, 111);
1143 // if ( mymat17->Determinant() == 0. ){
1144 // printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1145 // } else {
1146 // Double_t det = 0.;
1147 // TMatrixD invmatrix = (TMatrixD)(mymat17->Invert(&det));
1148 // printf(" Mymat3 determinant is %f \n",det);
1149 // cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1117);
1150 // };
1151 // cf->WriteFullMatrix(mymat17, 117);
1152
1153
1154 //
1155 // cf->UnLoadFullMatrix(i);
1156 // cf->UnLoadFullNMatrix(i);
1157 // delete fmatrix;
1158 // delete fnmat;
1159 //
1160 };
1161 };
1162 //
1163 printf(" done, closing file and exiting\n");
1164 //
1165 };
1166 //
1167 cf->CloseMatrixFile();
1168 //
1169 cf->Delete();
1170 //
1171 }

  ViewVC Help
Powered by ViewVC 1.1.23