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