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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.8 - (hide annotations) (download)
Fri Jan 25 15:09:07 2008 UTC (17 years, 9 months ago) by mocchiut
Branch: MAIN
Changes since 1.7: +588 -180 lines 
it works
1 mocchiut 1.1 /**
2     * \file CaloFranzini.cpp
3     * \author Emiliano Mocchiutti (2007/07/18)
4     */
5     //
6     // headers
7     //
8     #include <CaloFranzini.h>
9     //--------------------------------------
10     /**
11     * Default constructor
12     */
13     CaloFranzini::CaloFranzini(){
14     Clear();
15     }
16    
17     CaloFranzini::CaloFranzini(PamLevel2 *l2p){
18     //
19 mocchiut 1.8 lfile = NULL;
20     ffile = NULL;
21 mocchiut 1.1 brig = NULL;
22 mocchiut 1.4 brigm = NULL;
23 mocchiut 1.1 nbin = 0;
24     //
25     L2 = l2p;
26     //
27     if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
28     //
29     // Default variables
30     //
31     debug = false;
32     dolong = true;
33     dofull = false;
34 mocchiut 1.2 sntr = 0;
35 mocchiut 1.5 OBT = 0;
36     PKT = 0;
37     atime = 0;
38 mocchiut 1.1 //
39 mocchiut 1.5 crig = false;
40 mocchiut 1.3 sel = true;
41     cont = false;
42     N = 0;
43     NC = 43;
44     //
45     mask18b = -1;
46     //
47 mocchiut 1.1 Clear();
48     //
49     }
50    
51     void CaloFranzini::Clear(){
52     //
53 mocchiut 1.5 longtzeta = 0.;
54     fulltzeta = 0.;
55 mocchiut 1.2 degfre = 0;
56 mocchiut 1.8 fdegfre = 0;
57 mocchiut 1.1 memset(estrip, 0, 2*22*96*sizeof(Float_t));
58 mocchiut 1.3 memset(qplane, 0, 43*sizeof(Float_t));
59 mocchiut 1.1 //
60 mocchiut 1.8 numneg = 0;
61     numpos = 0;
62     negfulltzeta = 0.;
63     posfulltzeta = 0.;
64     aveposvar = 0.;
65     avenegvar = 0.;
66     minsvalue = numeric_limits<UInt_t>::max();
67     maxsvalue = numeric_limits<UInt_t>::min();
68     //
69 mocchiut 1.1 }
70    
71     void CaloFranzini::Print(){
72     //
73     Process();
74     //
75     printf("========================================================================\n");
76     printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
77     printf(" debug [debug flag]:.. %i\n",debug);
78 mocchiut 1.8 printf(" long degree of freedom :.. %i\n",this->GetLongDegreeOfFreedom());
79 mocchiut 1.1 printf(" longtzeta :.. %f\n",longtzeta);
80 mocchiut 1.3 printf(" longtzeta normalized :.. %f\n",this->GetNormLongTZeta());
81 mocchiut 1.8 printf(" full degree of freedom :.. %i\n",this->GetFullDegreeOfFreedom());
82     printf(" fulltzeta :.. %f\n",fulltzeta);
83     printf(" fulltzeta normalized :.. %f\n",this->GetNormFullTZeta());
84     printf(" fulltzeta negative contribute :.. %f\n",negfulltzeta);
85     printf(" fulltzeta positive contribute :.. %f\n",posfulltzeta);
86     printf(" fulltzeta number of negatives :.. %i\n",numneg);
87     printf(" fulltzeta number of positives :.. %i\n",numpos);
88     printf(" fulltzeta minimum variance :.. %f\n",minsvalue);
89     printf(" fulltzeta maximum variance :.. %f\n",maxsvalue);
90     printf(" fulltzeta average positive var :.. %f\n",aveposvar);
91     printf(" fulltzeta average negative var :.. %f\n",avenegvar);
92 mocchiut 1.1 printf("========================================================================\n");
93     //
94     }
95    
96     void CaloFranzini::Delete(){
97 mocchiut 1.2 //
98 mocchiut 1.8 if ( ffile ) ffile->Close();
99     if ( lfile ) lfile->Close();
100 mocchiut 1.2 //
101 mocchiut 1.1 Clear();
102 mocchiut 1.2 //
103     }
104    
105 mocchiut 1.3 void CaloFranzini::SetNoWpreSampler(Int_t n){
106 mocchiut 1.5 Int_t nc2 = NC/2;
107     if ( NC >= 37 ) nc2 = (NC+1)/2;
108     if ( nc2+n < 23 ){
109 mocchiut 1.3 N = n;
110     } else {
111     printf(" ERROR! Calorimeter is made of 22 W planes\n");
112     printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
113     printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
114     NC = 43;
115     N = 0;
116     };
117     }
118    
119     void CaloFranzini::SetNoWcalo(Int_t n){
120     if ( N+n < 23 ){
121     NC = n*2;
122     if ( NC >37 ) NC--;
123     } else {
124     printf(" ERROR! Calorimeter is made of 22 W planes\n");
125     printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
126     printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
127     NC = 43;
128     N = 0;
129     };
130     }
131 mocchiut 1.2
132     void CaloFranzini::Process(){
133     this->Process(0);
134     }
135    
136     void CaloFranzini::Process(Int_t itr){
137     //
138     if ( !L2 ){
139     printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
140     printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
141     return;
142     };
143     //
144 mocchiut 1.8 if ( !nbin || !brig || (!lfile && !ffile) ){
145 mocchiut 1.2 printf(" ERROR: it seems covariance matrix file has not been opened (CaloFranzini::Open()) \n");
146     printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
147     return;
148     };
149     //
150     Bool_t newentry = false;
151     //
152     if ( L2->IsORB() ){
153     if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || itr != sntr ){
154     newentry = true;
155     OBT = L2->GetOrbitalInfo()->OBT;
156     PKT = L2->GetOrbitalInfo()->pkt_num;
157     atime = L2->GetOrbitalInfo()->absTime;
158     sntr = itr;
159     };
160     } else {
161     newentry = true;
162     };
163     //
164     if ( !newentry ) return;
165     //
166     // Some variables
167     //
168     if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
169     //
170     this->Clear();
171     //
172     Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
173 mocchiut 1.5 if ( crig ) rig = L2->GetCaloLevel2()->qtot/260.;
174 mocchiut 1.2 //
175     // Fill the estrip matrix
176     //
177     Int_t nplane = 0;
178     Int_t view = 0;
179     Int_t plane = 0;
180     Int_t strip = 0;
181     Float_t mip = 0.;
182     //
183     //
184     //
185     if ( dolong ){
186     //
187 mocchiut 1.8 for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
188     //
189     mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
190     //
191     estrip[view][plane][strip] = mip;
192     //
193     nplane = 1 - view + 2 * (plane - N);
194     if ( nplane > (37-(2*N)) ) nplane--;
195     //
196     // if ( plane == (18+N) ) mip = 0.;
197     if ( nplane > -1 ) qplane[nplane] += mip;
198     //
199     };
200     //
201 mocchiut 1.3 if ( cont ){
202     for (Int_t i=0; i<22; i++){
203     if ( i == (18+N) ){
204     mask18b = 1 + 2 * (i - N);
205     break;
206     };
207     };
208     };
209     //
210     if ( sel ){
211     for (Int_t i=0; i<22; i++){
212     if ( i == (18-N) ){
213     mask18b = 1 + 2 * (i - N);
214     break;
215     };
216     };
217     };
218     //
219     if ( mask18b == 37 ) mask18b = -1;
220     //
221     Int_t dgf = 43;
222 mocchiut 1.2 //
223     for (Int_t i=0; i < 22; i++){
224     if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
225 mocchiut 1.3 dgf = 2 * i;
226 mocchiut 1.2 break;
227     };
228     if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
229 mocchiut 1.3 dgf = 1 + 2 * i;
230 mocchiut 1.2 break;
231     };
232     };
233     //
234 mocchiut 1.3 if ( dgf < 43 && dgf > 37 ) dgf--;
235     //
236     degfre = TMath::Min(dgf,NC);
237     //
238 mocchiut 1.8 // Float_t longzdiag = 0.;
239     // Float_t longzout = 0.;
240 mocchiut 1.5 //
241 mocchiut 1.2 if ( degfre > 0 ){
242     for (Int_t i = 0; i < degfre; i++){
243 mocchiut 1.3 if ( i != mask18b ){
244     for (Int_t j = 0; j < degfre; j++){
245     if ( j != mask18b ){
246 mocchiut 1.8 // if ( i == j ){
247     // longzdiag += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
248     // if ( debug ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig));
249     // } else {
250     // longzout += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
251     // if ( debug && i == (j+1) ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig));
252     // };
253     //
254 mocchiut 1.4 longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
255 mocchiut 1.8 //
256 mocchiut 1.3 };
257     };
258 mocchiut 1.2 };
259     };
260 mocchiut 1.8 // if ( debug ) printf(" diagonal total %f out of diagonal total %f \n",longzdiag,longzout);
261 mocchiut 1.2 };
262     };
263     if ( dofull ){
264 mocchiut 1.8 // printf(" ERROR: NOT IMPLEMENTED YET\n");
265     // printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
266     //
267     // FULL CALORIMETER
268     //
269     CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
270     //
271     Int_t dgf = 0;
272     Int_t cs = 0;
273     Int_t cd = 0;
274     Int_t mstrip = 0;
275     //
276     Int_t maxnpl = 0;
277     Float_t mipv[43][11];
278     memset(mipv,0,43*11*sizeof(Float_t));
279     //
280     for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
281     //
282     mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
283     //
284     // nplane = 1 - view + 2 * plane;
285     // if ( nplane > 37 ) nplane--;
286     nplane = 1 - view + 2 * (plane - N);
287     if ( nplane > (37-(2*N)) ) nplane--;
288     //
289     cs = ct->tibar[plane][view] - 1;
290     //
291     if ( ct->tibar[plane][view] > -1 ){
292     //
293     cd = 95 - cs;
294     //
295     mstrip = cd + strip;
296     //
297     Int_t lstr = this->ConvertStrip(mstrip);
298     //
299     if ( nplane > maxnpl ) maxnpl = nplane;
300     if ( nplane > -1 ) mipv[nplane][lstr] += mip;
301     //
302     };
303     // mipv[nplane][lstr] += mip;
304     //
305     };
306     //
307     Float_t mip1 = 1.;
308     Int_t cs1;
309     Int_t cd1;
310     Float_t mip2 = 1.;
311     Int_t cs2;
312     Int_t cd2;
313     Int_t mi = -1;
314     Int_t mj = -1;
315     Int_t nn1 = 0;
316     Int_t pl1 = 0;
317     Int_t vi1 = 0;
318     Int_t nn2 = 0;
319     Int_t pl2 = 0;
320     Int_t vi2 = 0;
321     Int_t mstrip1min = 0;
322     Int_t mstrip1max = 0;
323     Int_t mstrip2min = 0;
324     Int_t mstrip2max = 0;
325     //
326     Int_t tonpl = TMath::Min(maxnpl,NC);
327     //
328     Int_t rbi = 0;
329     for (Int_t i = 0; i<nbin-1; i++){
330     if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
331     rbi = i;
332     break;
333     };
334     };
335     //
336     //
337     Int_t therigb = rbi;
338     //
339     if ( rig < brigm->At(2) ){
340     // therigb = 0;
341     therigb = 2;
342     };
343     // if ( rig >= brigm->At(nbin-5) ){
344     if ( rig >= brigm->At(nbin-2) ){
345     therigb = nbin - 3;
346     // therigb = nbin - 5;
347     };
348     Int_t mtherig = therigb;
349     if ( mtherig >= 13 ) mtherig = 12;
350     //
351     if ( debug ) printf(" rig %f brigm0 %f brigm n2 %f nbin %i \n",rig,brigm->At(0),brigm->At(nbin-2),nbin);
352     //
353     if ( cont ){
354     for (Int_t i=0; i<22; i++){
355     if ( i == (18+N) ){
356     mask18b = 1 + 2 * (i - N);
357     break;
358     };
359     };
360     };
361     //
362     if ( sel ){
363     for (Int_t i=0; i<22; i++){
364     if ( i == (18-N) ){
365     mask18b = 1 + 2 * (i - N);
366     break;
367     };
368     };
369     };
370     //
371     if ( mask18b == 37 ) mask18b = -1;
372     //
373     if ( debug ) printf("Start main loop \n");
374     //
375     for (Int_t nplane1=0; nplane1<tonpl; nplane1++){
376     //
377     if ( nplane1 != mask18b ){
378     //
379     if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
380     vi1 = 1;
381     if ( nn1%2 ) vi1 = 0;
382     pl1 = (nn1 - 1 + vi1)/2;
383     //
384     cs1 = ct->tibar[pl1][vi1] - 1;
385     //
386     if ( ct->tibar[pl1][vi1] > -1 ){
387     //
388     dgf++;
389     //
390     cd1 = 95 - cs1;
391     //
392     Int_t at1 = TMath::Max(0,(cd1+0));
393     Int_t at2 = TMath::Min(190,(cd1+95));
394     mstrip1min = this->ConvertStrip(at1);
395     mstrip1max = this->ConvertStrip(at2) + 1;
396     //
397     for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
398     //
399     mj = -1;
400     //
401     // if ( debug ) printf(" get mip1 %i %i %f %i \n",nplane1,mstrip1,rig,therigb);
402     mip1 = mipv[nplane1][mstrip1] - this->GetFullAverageAt(nplane1,mstrip1,rig,therigb);
403     if ( mip1 < 0 ){
404     numneg++;
405     avenegvar += mip1;
406     };
407     if ( mip1 >= 0 ){
408     numpos++;
409     aveposvar += mip1;
410     };
411     if ( mip1 < minsvalue ) minsvalue = mip1;
412     if ( mip1 >= maxsvalue ) maxsvalue = mip1;
413     //
414     mi = (nplane1 * 11) + mstrip1;
415     //
416     if ( mip1 != 0. ){
417     //
418     for (Int_t nplane2=0; nplane2<tonpl; nplane2++){
419     //
420     if ( nplane2 != mask18b ){
421     //
422     if ( nplane2 >= 37 ) nn2 = nplane2 + 1;
423     vi2 = 1;
424     if ( nn2%2 ) vi2 = 0;
425     pl1 = (nn2 - 1 + vi2)/2;
426     //
427     cs2 = ct->tibar[pl2][vi2] - 1;
428     //
429     if ( ct->tibar[pl2][vi2] > -1 ){
430     //
431     cd2 = 95 - cs2;
432     //
433     Int_t t1 = TMath::Max(0,(cd2+0));
434     Int_t t2 = TMath::Min(190,(cd2+95));
435     mstrip2min = this->ConvertStrip(t1);
436     mstrip2max = this->ConvertStrip(t2) + 1;
437     //
438     for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
439     //
440     mip2 = mipv[nplane2][mstrip2] - this->GetFullAverageAt(nplane2,mstrip2,rig,therigb);
441     //
442     if ( mip2 != 0. ){
443     //
444     // Int_t sh = -5 + nplane1;
445     // if ( sh > 5 ) sh -= 11*nplane1;
446     // //
447     // mj = (nplane2 * 11) + mstrip2 + sh;
448     // //
449     // if ( mj < 0 ) mj += 473;
450     // if ( mj >= 473 ) mj -= 473;
451     // //
452     mj = (nplane2 * 11) + mstrip2;
453     //
454     Float_t svalue = mip1 * this->GetFullHmatrixAt(mi,mj,rig,mtherig,therigb) * mip2;
455     // fulltzeta += mip1 * this->GetFullHmatrixAt(mi,mj,rig,therigb) * mip2;
456     fulltzeta += svalue;
457     if ( svalue < 0. ) negfulltzeta += svalue;
458     if ( svalue > 0. ) posfulltzeta += svalue;
459     // (*fmatrix[rbi])[mi][mj] += (mip1 * mip2);
460     // if ( mstrip1 == mstrip2 ) printf("\n\n=> nplane1 %i nplane2 %i mstrip1 %i mstrip2 %i \n => mip1 %f H %f mip2 %f \n => mipv1 %f ave1 %f mipv2 %f ave2 %f \n => rig %f rigbin %i mtherigb %i\n",nplane1,nplane2,mstrip1,mstrip2,mip1,this->GetFullHmatrixAt(mi,mj,rig,mtherig),mip2,mipv[nplane1][mstrip1],this->GetFullAverageAt(nplane1,mstrip1,rig,therigb),mipv[nplane2][mstrip2],this->GetFullAverageAt(nplane2,mstrip2,rig,therigb),rig,therigb,mtherig);
461     //
462     };
463     };
464     };
465     };
466     };
467     };
468     };
469     };
470     };
471     };
472     //
473     fdegfre = dgf*11;
474     if ( numpos ) aveposvar /= numpos;
475     if ( numneg ) avenegvar /= numneg;
476     //
477     };
478 mocchiut 1.1 //
479 mocchiut 1.8 // if ( debug ) this->Print();
480 mocchiut 1.1 if ( debug ) printf(" exit \n");
481 mocchiut 1.2 }
482    
483    
484     Float_t CaloFranzini::GetNormLongTZeta(){
485     Process();
486     Float_t normz = 0.;
487 mocchiut 1.8 if ( degfre != 0 ) normz = longtzeta/(Float_t)degfre;
488 mocchiut 1.2 return normz;
489     }
490    
491     Float_t CaloFranzini::GetNormFullTZeta(){
492     Process();
493     Float_t normz = 0.;
494 mocchiut 1.8 if ( fdegfre != 0 ) normz = fulltzeta/(Float_t)fdegfre;
495 mocchiut 1.2 return normz;
496 mocchiut 1.1 }
497    
498     Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){
499     //
500 mocchiut 1.8 lfile = new TFile(matrixfile.Data(),"READ");
501 mocchiut 1.1 //
502 mocchiut 1.8 if ( !lfile || lfile->IsZombie() ){
503     if ( dolong ){
504     lfile = new TFile(matrixfile.Data(),"RECREATE");
505     printf(" Create file %s \n",lfile->GetName());
506     };
507     if ( dofull ){
508     ffile = new TFile(matrixfile.Data(),"RECREATE");
509     printf(" Create file %s \n",ffile->GetName());
510     };
511 mocchiut 1.1 } else {
512 mocchiut 1.8 lfile->Close();
513 mocchiut 1.1 printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
514     return(false);
515     };
516     //
517     return(true);
518     //
519     }
520    
521 mocchiut 1.4 Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){
522     //
523 mocchiut 1.8 if ( dolong ){
524     lfile = new TFile(matrixfile.Data(),"UPDATE");
525     //
526     if ( !lfile || lfile->IsZombie() ){
527     printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
528     return(false);
529     };
530     };
531 mocchiut 1.4 //
532 mocchiut 1.8 if ( dofull ){
533     ffile = new TFile(matrixfile.Data(),"UPDATE");
534     //
535     if ( !ffile || ffile->IsZombie() ){
536     printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
537     return(false);
538     };
539 mocchiut 1.4 };
540     //
541     return(true);
542     //
543     }
544    
545 mocchiut 1.1 void CaloFranzini::WriteNumBin(Int_t numbin){
546 mocchiut 1.8 if ( dolong ){
547     lfile->cd();
548     TArrayI nbi(1, &numbin);
549     lfile->WriteObject(&nbi, "nbinenergy");
550     };
551     if ( dofull ){
552     ffile->cd();
553     TArrayI nbi(1, &numbin);
554     ffile->WriteObject(&nbi, "nbinenergy");
555     };
556 mocchiut 1.1 }
557    
558     void CaloFranzini::WriteRigBin(TArrayF *rigbin){
559 mocchiut 1.8 if ( dolong ){
560     lfile->cd();
561     // rigbin->Write("binrig");
562     lfile->WriteObject(&(*rigbin), "binrig");
563     };
564     if ( dofull ){
565     ffile->cd();
566     // rigbin->Write("binrig");
567     ffile->WriteObject(&(*rigbin), "binrig");
568     };
569 mocchiut 1.1 }
570    
571     void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
572 mocchiut 1.8 lfile->cd();
573 mocchiut 1.1 TString name = Form("qplmeann%i",bin);
574 mocchiut 1.8 lfile->WriteObject(&(*qpl),name.Data());
575 mocchiut 1.3 }
576    
577 mocchiut 1.5 void CaloFranzini::WriteFullMean(TMatrixD *qpl, Int_t bin){
578 mocchiut 1.8 ffile->cd();
579 mocchiut 1.5 TString name = Form("fqplmeann%i",bin);
580 mocchiut 1.8 ffile->WriteObject(&(*qpl),name.Data());
581     ffile->Purge();
582 mocchiut 1.6 }
583    
584     void CaloFranzini::WriteFullNMean(TMatrixD *qpl, Int_t bin){
585 mocchiut 1.8 ffile->cd();
586 mocchiut 1.6 TString name = Form("fnqplmeann%i",bin);
587 mocchiut 1.8 ffile->WriteObject(&(*qpl),name.Data());
588     ffile->Purge();
589 mocchiut 1.5 }
590    
591 mocchiut 1.3 void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
592 mocchiut 1.8 lfile->cd();
593 mocchiut 1.3 TString name = Form("matrixn%i",bin);
594     // mat.Write(name.Data());
595 mocchiut 1.8 lfile->WriteObject(&mat,name.Data());
596 mocchiut 1.1 }
597    
598 mocchiut 1.7 void CaloFranzini::WriteInvertedFullMatrix(TMatrixD mat, Int_t bin){
599 mocchiut 1.8 ffile->cd();
600 mocchiut 1.5 TString name = Form("fmatrixn%i",bin);
601     // mat.Write(name.Data());
602 mocchiut 1.8 ffile->WriteObject(&mat,name.Data());
603 mocchiut 1.5 }
604    
605 mocchiut 1.1 void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
606 mocchiut 1.8 lfile->cd();
607 mocchiut 1.3 TString name = Form("origmatrixn%i",bin);
608     // mat.Write(name.Data());
609 mocchiut 1.8 lfile->WriteObject(&(*mat),name.Data());
610 mocchiut 1.1 }
611    
612 mocchiut 1.7 void CaloFranzini::WriteFullMatrix(TMatrixD *mat, Int_t bin){
613 mocchiut 1.8 ffile->cd();
614 mocchiut 1.5 TString name = Form("origfmatrixn%i",bin);
615     // mat.Write(name.Data());
616 mocchiut 1.8 ffile->WriteObject(&(*mat),name.Data());
617     ffile->Purge();
618 mocchiut 1.6 }
619    
620     void CaloFranzini::WriteFullNMatrix(TMatrixF *mat, Int_t bin){
621 mocchiut 1.8 ffile->cd();
622 mocchiut 1.6 TString name = Form("origfnmatrixn%i",bin);
623     // mat.Write(name.Data());
624 mocchiut 1.8 ffile->WriteObject(&(*mat),name.Data());
625     ffile->Purge();
626 mocchiut 1.5 }
627    
628 mocchiut 1.1 void CaloFranzini::CloseMatrixFile(){
629 mocchiut 1.8 if ( dolong ){
630     lfile->cd();
631     lfile->Close();
632     };
633     if ( dofull ){
634     ffile->cd();
635     ffile->Close();
636     };
637 mocchiut 1.1 }
638    
639 mocchiut 1.8 Bool_t CaloFranzini::Open(TString matrixfile){
640     return(this->Open(matrixfile,""));
641     }
642 mocchiut 1.1
643 mocchiut 1.8 Bool_t CaloFranzini::Open(TString longmatrixfile, TString fullmatrixfile){
644 mocchiut 1.1 //
645     // find matrix file
646     //
647 mocchiut 1.8 if ( !strcmp(longmatrixfile.Data(),"") ){
648     if (dolong) longmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_longel.root";
649 mocchiut 1.1 };
650 mocchiut 1.8 if ( !strcmp(fullmatrixfile.Data(),"") ){
651     if (dofull) fullmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";
652 mocchiut 1.4 };
653 mocchiut 1.5 //
654     if ( dolong ){
655 mocchiut 1.8 //
656     lfile = new TFile(longmatrixfile.Data(),"READ");
657     //
658     if ( !lfile || lfile->IsZombie() ){
659     printf(" ERROR: cannot open file %s \n",longmatrixfile.Data());
660     return(false);
661     };
662     //
663     if ( !this->LoadBin() ){
664     printf(" %s \n",longmatrixfile.Data());
665     return(false);
666     };
667     //
668 mocchiut 1.5 if ( !this->LoadLong() ){
669 mocchiut 1.8 printf(" %s \n",longmatrixfile.Data());
670 mocchiut 1.5 return(false);
671     };
672     //
673     if ( !this->LoadMatrices() ){
674 mocchiut 1.8 printf(" %s \n",longmatrixfile.Data());
675 mocchiut 1.5 return(false);
676     };
677     };
678     //
679     if ( dofull ){
680 mocchiut 1.8 //
681     ffile = new TFile(fullmatrixfile.Data(),"READ");
682     //
683     if ( !ffile || ffile->IsZombie() ){
684     printf(" ERROR: cannot open file %s \n",fullmatrixfile.Data());
685     return(false);
686     };
687     //
688     if ( !dolong ){
689     //
690     if ( !this->LoadBin(true) ){
691     printf(" %s \n",fullmatrixfile.Data());
692     return(false);
693     };
694     //
695     };
696 mocchiut 1.5 if ( !this->LoadFull() ){
697 mocchiut 1.8 printf(" %s \n",fullmatrixfile.Data());
698 mocchiut 1.5 return(false);
699     };
700     //
701     if ( !this->LoadFullMatrices() ){
702 mocchiut 1.8 printf(" %s \n",fullmatrixfile.Data());
703 mocchiut 1.5 return(false);
704     };
705 mocchiut 1.4 };
706     //
707     //
708     return(true);
709     //
710     }
711    
712 mocchiut 1.8
713 mocchiut 1.4 Bool_t CaloFranzini::LoadBin(){
714 mocchiut 1.8 return(this->LoadBin(false));
715     }
716    
717     Bool_t CaloFranzini::LoadBin(Bool_t full){
718 mocchiut 1.4 //
719 mocchiut 1.8 if ( !full ) {
720     //
721     TArrayI *numbin = (TArrayI*)lfile->Get("nbinenergy");
722     if ( !numbin ){
723     printf(" ERROR: cannot read number of bins from file ");
724     return(false);
725     };
726     nbin = (Int_t)numbin->At(0);
727     if ( nbin <= 0 ){
728     printf(" ERROR: cannot work with 0 energy bins from file ");
729     return(false);
730     };
731     //
732     brig = (TArrayF*)lfile->Get("binrig");
733     if ( !brig ){
734     printf(" ERROR: cannot read rigidity binning from file ");
735     return(false);
736     };
737     //
738     brigm=(TArrayF*)lfile->Get("binrigmean");
739     if ( !brigm ){
740     printf(" ERROR: cannot read mean rigidity binning from file ");
741     return(false);
742     };
743     //
744     } else {
745     //
746     TArrayI *numbin = (TArrayI*)ffile->Get("nbinenergy");
747     if ( !numbin ){
748     printf(" ERROR: cannot read number of bins from file ");
749     return(false);
750     };
751     nbin = (Int_t)numbin->At(0);
752     if ( nbin <= 0 ){
753     printf(" ERROR: cannot work with 0 energy bins from file ");
754     return(false);
755     };
756     //
757     brig = (TArrayF*)ffile->Get("binrig");
758     if ( !brig ){
759     printf(" ERROR: cannot read rigidity binning from file ");
760     return(false);
761     };
762     //
763     brigm=(TArrayF*)ffile->Get("binrigmean");
764     if ( !brigm ){
765     printf(" ERROR: cannot read mean rigidity binning from file ");
766     return(false);
767     };
768 mocchiut 1.4 };
769     //
770 mocchiut 1.5 return(true);
771     }
772    
773     Bool_t CaloFranzini::LoadLong(){
774     //
775 mocchiut 1.4 for (Int_t i=0;i<17;i++){
776 mocchiut 1.5 TString name = Form("qplmeann%i",i);
777 mocchiut 1.8 qplmean[i] = (TArrayF*)lfile->Get(name.Data());
778 mocchiut 1.5 if ( !qplmean[i] ){
779     printf(" ERROR: cannot read average from file ");
780     return(false);
781     };
782     };
783     //
784     return(true);
785     }
786    
787     Bool_t CaloFranzini::LoadFull(){
788     //
789     for (Int_t i=0;i<17;i++){
790     TString name = Form("fqplmeann%i",i);
791 mocchiut 1.8 fqplmean[i] = (TMatrixD*)ffile->Get(name.Data());
792 mocchiut 1.5 if ( !fqplmean[i] ){
793     printf(" ERROR: cannot read average from file ");
794     return(false);
795     };
796 mocchiut 1.4 };
797     //
798 mocchiut 1.1 return(true);
799 mocchiut 1.4 }
800    
801     Bool_t CaloFranzini::LoadMatrices(){
802 mocchiut 1.1 //
803 mocchiut 1.4 for (Int_t i=0;i<17;i++){
804     TString name1 = Form("matrixn%i",i);
805 mocchiut 1.8 hmat[i] = (TMatrixD*)lfile->Get(name1.Data());
806 mocchiut 1.4 };
807     //
808     return(true);
809 mocchiut 1.1 }
810    
811 mocchiut 1.5 Bool_t CaloFranzini::LoadFullMatrices(){
812     //
813     for (Int_t i=0;i<17;i++){
814     TString name1 = Form("fmatrixn%i",i);
815 mocchiut 1.8 hfmat[i] = (TMatrixD*)ffile->Get(name1.Data());
816 mocchiut 1.5 };
817     //
818     return(true);
819     }
820    
821 mocchiut 1.1 TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
822     //
823 mocchiut 1.4 Int_t mv = 0;
824 mocchiut 1.1 for (Int_t i = 0; i<nbin-1; i++){
825     if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
826 mocchiut 1.4 mv = i;
827 mocchiut 1.1 break;
828     };
829     };
830     if ( rig < brig->At(0) ){
831     printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
832 mocchiut 1.4 mv = 0;
833 mocchiut 1.1 };
834 mocchiut 1.3 if ( rig >= brig->At(nbin-1) ){
835     printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
836 mocchiut 1.4 mv = nbin-2;
837 mocchiut 1.1 };
838     //
839 mocchiut 1.4 return(hmat[mv]);
840 mocchiut 1.1 //
841     }
842    
843 mocchiut 1.6 TMatrixD *CaloFranzini::LoadFullAverage(Int_t rigbin){
844     //
845     TString name = Form("fqplmeann%i",rigbin);
846 mocchiut 1.8 TMatrixD *fmean=(TMatrixD*)ffile->Get(name.Data());
847 mocchiut 1.6 //
848     return(fmean);
849     //
850     }
851    
852     TMatrixF *CaloFranzini::LoadFullMatrix(Int_t rigbin){
853     //
854     TString name = Form("origfmatrixn%i",rigbin);
855 mocchiut 1.8 TMatrixF *fmatri=(TMatrixF*)ffile->Get(name.Data());
856 mocchiut 1.6 //
857     return(fmatri);
858     //
859     }
860    
861     void CaloFranzini::LoadFullMatrix(Int_t rigbin, TMatrixF *&fmatri){
862     //
863     TString name = Form("origfmatrixn%i",rigbin);
864 mocchiut 1.8 fmatri=(TMatrixF*)ffile->Get(name.Data());
865 mocchiut 1.6 //
866     }
867    
868     void CaloFranzini::UnLoadFullAverage(Int_t rigbin){
869     //
870     TString name = Form("fqplmeann%i",rigbin);
871 mocchiut 1.8 ffile->Delete(name.Data());
872 mocchiut 1.6 //
873     }
874    
875     void CaloFranzini::UnLoadFullMatrix(Int_t rigbin){
876     //
877     TString name = Form("origfmatrixn%i",rigbin);
878 mocchiut 1.8 ffile->Delete(name.Data());
879 mocchiut 1.6 //
880     }
881    
882     TMatrixF *CaloFranzini::LoadFullNMatrix(Int_t rigbin){
883     //
884     TString name = Form("origfnmatrixn%i",rigbin);
885 mocchiut 1.8 TMatrixF *fnmatri=(TMatrixF*)ffile->Get(name.Data());
886 mocchiut 1.6 //
887     return(fnmatri);
888     //
889     }
890    
891     void CaloFranzini::UnLoadFullNMatrix(Int_t rigbin){
892     //
893     TString name = Form("origfnmatrixn%i",rigbin);
894 mocchiut 1.8 ffile->Delete(name.Data());
895 mocchiut 1.6 //
896     }
897    
898     TMatrixD *CaloFranzini::LoadFullNAverage(Int_t rigbin){
899     //
900     TString name = Form("fnqplmeann%i",rigbin);
901 mocchiut 1.8 TMatrixD *fnmean=(TMatrixD*)ffile->Get(name.Data());
902 mocchiut 1.6 //
903     return(fnmean);
904     //
905     }
906    
907     void CaloFranzini::UnLoadFullNAverage(Int_t rigbin){
908     //
909     TString name = Form("fnqplmeann%i",rigbin);
910 mocchiut 1.8 ffile->Delete(name.Data());
911 mocchiut 1.6 //
912     }
913    
914 mocchiut 1.1
915     TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
916     //
917 mocchiut 1.4 Int_t mv=0;
918 mocchiut 1.1 for (Int_t i = 0; i<nbin-1; i++){
919     if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
920 mocchiut 1.4 mv = i;
921 mocchiut 1.1 break;
922     };
923     };
924     if ( rig < brig->At(0) ){
925     printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
926 mocchiut 1.4 mv = 0;
927 mocchiut 1.1 };
928 mocchiut 1.3 if ( rig >= brig->At(nbin-1) ){
929     printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
930 mocchiut 1.4 mv=nbin-2;
931 mocchiut 1.1 };
932     //
933 mocchiut 1.4 return(qplmean[mv]);
934 mocchiut 1.1 //
935 mocchiut 1.4 }
936    
937     Float_t CaloFranzini::GetAverageAt(Int_t plane, Float_t rig){
938     //
939     Int_t therigb = 0;
940     for (Int_t i = 0; i<nbin-2; i++){
941     if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
942     therigb = i;
943     break;
944     };
945     };
946     //
947     Float_t minrig;
948     Float_t maxrig;
949     //
950     //
951     if ( rig < brigm->At(0) ){
952     if ( rig < brig->At(0) ){
953 mocchiut 1.5 // printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
954 mocchiut 1.4 };
955     therigb = 0;
956     };
957     if ( rig >= brigm->At(nbin-2) ){
958     if ( rig >= brig->At(nbin-2) ) {
959 mocchiut 1.5 // printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
960 mocchiut 1.4 };
961     therigb = nbin - 3;
962     };
963     //
964     minrig = brigm->At(therigb);
965     maxrig = brigm->At(therigb+1);
966     //
967     Float_t minene = (*qplmean[therigb])[plane];
968     Float_t maxene = (*qplmean[therigb+1])[plane];
969     //
970     if ( maxrig == minrig ){
971     printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
972     return(0.);
973     };
974     Float_t b = log(maxene/minene)/(maxrig-minrig);
975     Float_t a = minene/exp(minrig*b);
976     Float_t ave = a*exp(b*rig);
977 mocchiut 1.6 if ( a == 0. || minene == 0. || ave != ave ){
978     // if ( a == 0. || minene == 0. ){
979     Float_t m = (maxene-minene)/(maxrig-minrig);
980     Float_t q = minene - m * minrig;
981     ave = rig * m + q;
982     };
983 mocchiut 1.4 //
984     return(ave);
985 mocchiut 1.1 //
986     }
987 mocchiut 1.6
988 mocchiut 1.5 Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig){
989     //
990     Int_t therigb = 0;
991     for (Int_t i = 0; i<nbin-2; i++){
992     if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
993     therigb = i;
994     break;
995     };
996     };
997     //
998     //
999     if ( rig < brigm->At(0) ){
1000     if ( rig < brig->At(0) ){
1001     // printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
1002     };
1003     therigb = 0;
1004     };
1005     if ( rig >= brigm->At(nbin-2) ){
1006     if ( rig >= brig->At(nbin-2) ) {
1007     // printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
1008     };
1009     therigb = nbin - 3;
1010     };
1011     //
1012 mocchiut 1.6 return(this->GetFullAverageAt(plane,strip,rig,therigb));
1013     //
1014     }
1015    
1016     Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig, Int_t therigb){
1017     //
1018     Float_t minrig;
1019     Float_t maxrig;
1020     //
1021 mocchiut 1.5 minrig = brigm->At(therigb);
1022     maxrig = brigm->At(therigb+1);
1023     //
1024     Float_t minene = (*fqplmean[therigb])[plane][strip];
1025     Float_t maxene = (*fqplmean[therigb+1])[plane][strip];
1026     //
1027     if ( maxrig == minrig ){
1028     printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
1029     return(0.);
1030     };
1031     Float_t b = log(maxene/minene)/(maxrig-minrig);
1032     Float_t a = minene/exp(minrig*b);
1033     Float_t ave = a*exp(b*rig);
1034 mocchiut 1.6 if ( a == 0. || minene == 0. || ave != ave ){
1035     Float_t m = (maxene-minene)/(maxrig-minrig);
1036     Float_t q = minene - m * minrig;
1037 mocchiut 1.8 ave = rig * m + q;
1038 mocchiut 1.6 };
1039     //
1040     // ave += (44.-plane)*strip;
1041     //if ( a == 0. ) ave = 0.;
1042     if ( !(ave == ave) ) printf("a %f b %f ave %f maxene %f minene %f maxrig %f minrig %f \n",a,b,ave,maxene,minene,maxrig,minrig);
1043 mocchiut 1.5 //
1044     return(ave);
1045     //
1046     }
1047 mocchiut 1.3
1048 mocchiut 1.4 Float_t CaloFranzini::GetHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
1049     Int_t therigb = 0;
1050     for (Int_t i = 0; i<nbin-2; i++){
1051     if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
1052     therigb = i;
1053     break;
1054     };
1055     };
1056     //
1057     Float_t minrig;
1058     Float_t maxrig;
1059     //
1060     if ( rig < brigm->At(0) ){
1061     if ( rig < brig->At(0) ){
1062 mocchiut 1.8 // printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
1063 mocchiut 1.4 };
1064     therigb = 0;
1065     };
1066 mocchiut 1.8 // if ( rig >= brigm->At(nbin-4) ){ // -2
1067 mocchiut 1.5 if ( rig >= brigm->At(nbin-2) ){
1068 mocchiut 1.4 if ( rig >= brig->At(nbin-2) ) {
1069 mocchiut 1.8 // printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
1070 mocchiut 1.4 };
1071 mocchiut 1.8 // therigb = nbin - 5;// -3
1072 mocchiut 1.4 therigb = nbin - 3;
1073     };
1074     //
1075 mocchiut 1.5 if ( therigb < 2 ) therigb = 2;
1076 mocchiut 1.4 minrig = brigm->At(therigb);
1077     maxrig = brigm->At(therigb+1);
1078 mocchiut 1.5 // printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
1079 mocchiut 1.4 //
1080     Float_t minene = (*hmat[therigb])[iindex][jindex];
1081     Float_t maxene = (*hmat[therigb+1])[iindex][jindex];
1082 mocchiut 1.5 // printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,a,b,rig,ave);
1083     //
1084     // Float_t a = 0.;
1085     // Float_t b = 0.;
1086     // Float_t ave = 0.;
1087     // if ( minene == 0. ){
1088     //
1089     // } else {
1090     // b = log(maxene/minene)/(maxrig-minrig);
1091     // a = minene/exp(minrig*b);
1092     // ave = a*exp(b*rig);
1093     // };
1094     //
1095     Float_t m = (maxene-minene)/(maxrig-minrig);
1096     Float_t q = minene - m * minrig;
1097     Float_t ave = rig * m + q;
1098     if ( debug ) printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,m,q,rig,ave);
1099     //
1100     //
1101     return(ave);
1102     //
1103     }
1104    
1105     Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
1106     Int_t therigb = 0;
1107     for (Int_t i = 0; i<nbin-2; i++){
1108     if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
1109     therigb = i;
1110     break;
1111     };
1112     };
1113     //
1114     if ( rig < brigm->At(0) ){
1115     if ( rig < brig->At(0) ){
1116 mocchiut 1.8 // printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
1117 mocchiut 1.5 };
1118     therigb = 0;
1119     };
1120     // if ( rig >= brigm->At(nbin-4) ){ // -2
1121     if ( rig >= brigm->At(nbin-2) ){
1122 mocchiut 1.8 //if ( rig >= brigm->At(nbin-5) ){
1123 mocchiut 1.5 if ( rig >= brig->At(nbin-2) ) {
1124 mocchiut 1.8 // printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
1125 mocchiut 1.5 };
1126 mocchiut 1.8 // therigb = nbin - 5;// -3
1127     // therigb = nbin - 5;// -3
1128 mocchiut 1.5 therigb = nbin - 3;
1129     };
1130     //
1131     if ( therigb < 2 ) therigb = 2;
1132 mocchiut 1.8 if ( therigb > 13 ) therigb = 13;
1133 mocchiut 1.6 //
1134     return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb));
1135     //
1136     }
1137    
1138 mocchiut 1.8
1139 mocchiut 1.6 Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb){
1140     //
1141 mocchiut 1.8 return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb,0));
1142     //
1143     }
1144    
1145     Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb, Int_t mtherig){
1146     //
1147     // Int_t lofit = 12;
1148     Int_t lofit = 3;
1149     //
1150     Float_t lowrig = 0.;
1151 mocchiut 1.6 Float_t minrig;
1152     Float_t maxrig;
1153     //
1154 mocchiut 1.8 if ( mtherig > lofit ) lowrig = brigm->At(therigb-1);
1155 mocchiut 1.5 minrig = brigm->At(therigb);
1156     maxrig = brigm->At(therigb+1);
1157 mocchiut 1.8 //if ( therigb > 10 ) printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
1158 mocchiut 1.5 //
1159 mocchiut 1.8 Float_t lowene = 0.;
1160     if ( mtherig > lofit ) lowene = (*hfmat[therigb-1])[iindex][jindex];
1161 mocchiut 1.5 Float_t minene = (*hfmat[therigb])[iindex][jindex];
1162     Float_t maxene = (*hfmat[therigb+1])[iindex][jindex];
1163     // printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,a,b,rig,ave);
1164     //
1165 mocchiut 1.8 Float_t ave = 0.;
1166     //
1167     // Float_t a = 0.;
1168     // Float_t b = 0.;
1169     // Float_t ave = 0.;
1170     // if ( minene == 0. ){
1171     //
1172     // } else {
1173     // b = log(maxene/minene)/(maxrig-minrig);
1174     // a = minene/exp(minrig*b);
1175     // ave = a*exp(b*rig);
1176     // };
1177 mocchiut 1.5 //
1178 mocchiut 1.8 if ( mtherig > lofit ){
1179     //
1180     // FIT
1181     //
1182     // Float_t x[3]={lowrig,minrig,maxrig};
1183     // Float_t y[3]={lowene,minene,maxene};
1184     // //
1185     // TGraph *tg= new TGraph(3,x,y);
1186     // //
1187     // // gStyle->SetLabelSize(0.04);
1188     // // gStyle->SetNdivisions(510,"XY");
1189     // // TCanvas *c = new TCanvas();
1190     // // c->Draw();
1191     // // c->cd();
1192     // // tg->Draw("AP");
1193     // //
1194     // TF1 *fun = new TF1("fun","pol2");
1195     // tg->Fit("fun","QNC");
1196     // //
1197     // ave = fun->GetParameter(0) + rig * fun->GetParameter(1) + rig * rig * fun->GetParameter(2);
1198     // //
1199     // // printf(" therigb %i ave %f rig %f lowrig %f minrig %f maxrig %f lowene %f minene %f maxene %f \n",therigb,ave,rig,lowrig,minrig,maxrig,lowene,minene,maxene);
1200     // //
1201     // //
1202     // // c->Modified();
1203     // // c->Update();
1204     // // gSystem->ProcessEvents();
1205     // // gSystem->Sleep(6000);
1206     // // c->Close();
1207     // // gStyle->SetLabelSize(0);
1208     // // gStyle->SetNdivisions(1,"XY");
1209     // //
1210     // tg->Delete();
1211     // fun->Delete();
1212     Float_t mrigl = (lowrig+minrig)/2.;
1213     Float_t mrigu = (minrig+maxrig)/2.;
1214     Float_t menel = (lowene+minene)/2.;
1215     Float_t meneu = (minene+maxene)/2.;
1216     Float_t m = (meneu-menel)/(mrigu-mrigl);
1217     Float_t q = menel - m * mrigl;
1218     ave = rig * m + q;
1219     //
1220     } else {
1221     Float_t m = (maxene-minene)/(maxrig-minrig);
1222     Float_t q = minene - m * minrig;
1223     ave = rig * m + q;
1224     if ( debug ) printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,m,q,rig,ave);
1225     };
1226 mocchiut 1.4 //
1227     //
1228     return(ave);
1229     //
1230     }
1231 mocchiut 1.3
1232     void CaloFranzini::DrawLongAverage(Float_t rig){
1233     //
1234     TArrayF *ll = this->LoadLongAverage(rig);
1235     //
1236     gStyle->SetLabelSize(0.04);
1237     gStyle->SetNdivisions(510,"XY");
1238     //
1239     TString hid = Form("cfralongvyvx");
1240     TCanvas *tcf = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1241     if ( tcf ){
1242     tcf->cd();
1243     } else {
1244     tcf = new TCanvas(hid,hid);
1245     };
1246     //
1247     TString thid = Form("hfralongvyvx");
1248     TH1F *thf = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1249     if ( thf ) thf->Delete();
1250     thf = new TH1F(thid,thid,44,-0.5,43.5);
1251     tcf->cd();
1252     Float_t qpl[43];
1253     for (Int_t st=0;st<43;st++){
1254     qpl[st] = ll->At(st);
1255 mocchiut 1.4 printf("st %i qpl %f\n",st,qpl[st]);
1256 mocchiut 1.3 };
1257     for (Int_t st=0;st<44;st++){
1258     if ( st == 37 ){
1259     thf->Fill(st,0.);
1260     } else {
1261     Int_t ss = st;
1262     if ( st > 37 ) ss--;
1263     thf->Fill(st,qpl[ss]);
1264     };
1265     };
1266     thf->Draw();
1267     tcf->Modified();
1268     tcf->Update();
1269     //
1270     gStyle->SetLabelSize(0);
1271     gStyle->SetNdivisions(1,"XY");
1272     //
1273     };
1274 mocchiut 1.5
1275     void CaloFranzini::DrawLongAverage(Int_t bin){
1276     //
1277     TArrayF *ll = this->LoadLongAverage(brigm->At(bin));
1278     //
1279     gStyle->SetLabelSize(0.04);
1280     gStyle->SetNdivisions(510,"XY");
1281     //
1282     TString hid = Form("cfralongvyvx");
1283     TCanvas *tcf = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1284     if ( tcf ){
1285     tcf->cd();
1286     } else {
1287     tcf = new TCanvas(hid,hid);
1288     };
1289     //
1290     TString thid = Form("hfralongvyvx");
1291     TH1F *thf = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1292     if ( thf ) thf->Delete();
1293     thf = new TH1F(thid,thid,44,-0.5,43.5);
1294     tcf->cd();
1295     Float_t qpl[43];
1296     for (Int_t st=0;st<43;st++){
1297     qpl[st] = ll->At(st);
1298     printf("st %i qpl %f\n",st,qpl[st]);
1299     };
1300     for (Int_t st=0;st<44;st++){
1301     if ( st == 37 ){
1302     thf->Fill(st,0.);
1303     } else {
1304     Int_t ss = st;
1305     if ( st > 37 ) ss--;
1306     thf->Fill(st,qpl[ss]);
1307     };
1308     };
1309     thf->Draw();
1310     tcf->Modified();
1311     tcf->Update();
1312     //
1313     gStyle->SetLabelSize(0);
1314     gStyle->SetNdivisions(1,"XY");
1315     //
1316     };
1317 mocchiut 1.7
1318     Int_t CaloFranzini::ConvertStrip(Int_t mstrip){
1319     //
1320     Int_t lastrip = 0;
1321 mocchiut 1.8 // //
1322     // if ( mstrip < 50 ) lastrip = 0;
1323     // //
1324     // if ( mstrip >= 50 && mstrip < 64 ) lastrip = 1;
1325     // //
1326     // if ( mstrip >= 64 && mstrip < 70 ) lastrip = 2;
1327     // //
1328     // if ( mstrip >= 70 && mstrip < 75 ) lastrip = 3;
1329     // //
1330     // if ( mstrip >= 75 && mstrip < 84 ){
1331     // lastrip = (int)trunc((mstrip - 75)/3.) + 4;
1332     // };
1333     // //
1334     // if ( mstrip >= 84 && mstrip < 90 ){
1335     // lastrip = (int)trunc((mstrip - 84)/2.) + 7;
1336     // };
1337     // //
1338     // if ( mstrip >= 90 && mstrip < 101 ){
1339     // lastrip = mstrip - 90 + 10;
1340     // };
1341     // //
1342     // if ( mstrip >= 101 && mstrip < 107 ){
1343     // lastrip = (int)trunc((mstrip - 101)/2.) + 21;
1344     // };
1345     // //
1346     // //
1347     // if ( mstrip >= 107 && mstrip < 116 ){
1348     // lastrip = (int)trunc((mstrip - 107)/3.) + 24;
1349     // };
1350     // //
1351     // if ( mstrip >= 116 && mstrip < 121 ) lastrip = 27;
1352     // //
1353     // if ( mstrip >= 121 && mstrip < 127 ) lastrip = 28;
1354     // //
1355     // if ( mstrip >= 127 && mstrip < 141 ) lastrip = 29;
1356     // //
1357     // if ( mstrip >= 141 ) lastrip = 30;
1358     // //
1359     //
1360     if ( mstrip < 83 ) lastrip = 0;
1361     //
1362     if ( mstrip >= 83 && mstrip < 90 ) lastrip = 1;
1363     //
1364     if ( mstrip >= 90 && mstrip < 93 ) lastrip = 2;
1365     //
1366     if ( mstrip >= 93 && mstrip < 98 ){
1367     lastrip = mstrip - 93 + 3;
1368 mocchiut 1.7 };
1369     //
1370 mocchiut 1.8 if ( mstrip >= 98 && mstrip < 101 ) lastrip = 8;
1371 mocchiut 1.7 //
1372 mocchiut 1.8 if ( mstrip >= 101 && mstrip < 107 ) lastrip = 9;
1373 mocchiut 1.7 //
1374 mocchiut 1.8 if ( mstrip >= 107 ) lastrip = 10;
1375 mocchiut 1.7 //
1376     return(lastrip);
1377     }
  ViewVC Help
Powered by ViewVC 1.1.23