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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by mocchiut, Thu Dec 13 17:08:09 2007 UTC revision 1.6 by mocchiut, Tue Jan 15 12:41:38 2008 UTC
# Line 6  Line 6 
6  #include <TH1F.h>  #include <TH1F.h>
7  #include <TH2F.h>  #include <TH2F.h>
8  #include <TMatrixD.h>  #include <TMatrixD.h>
9    #include <TMatrixF.h>
10  #include <TArrayF.h>  #include <TArrayF.h>
11  #include <TArrayI.h>  #include <TArrayI.h>
12  //  //
# Line 14  Line 15 
15  //  //
16  using namespace std;  using namespace std;
17  //  //
18    extern Bool_t MATRIX;
19    extern Bool_t FULL;
20    extern Bool_t SIMU;
21    extern Bool_t CRIG;
22    extern Bool_t SRIG;
23  CaloFranzini *cf;  CaloFranzini *cf;
24  Int_t nbin;  Int_t nbin;
25  Float_t rig[18];  Float_t rig[18];
26  Float_t rmean[18];  Float_t rmean[17];
27  Float_t qqplane[18][43];  Int_t ntot[17];
28  Int_t nnqplane[18][43];  Int_t MDIM = 8213;
29    //Int_t MDIM = 4128;
30  TArrayF *qplane[18];  //Float_t qqplane[17][43];
31  TArrayI *nqplane[18];  //Int_t nnqplane[17][43];
32  TMatrixD *matrix[18];  
33  TMatrixD *nmat[18];  TArrayF *qplane[17];
34    TArrayI *nqplane[17];
35    TMatrixD *matrix[17];
36    TMatrixD *nmat[17];
37    
38    TMatrixD *fqplane;
39    TMatrixD *fnqplane;
40    //TMatrixD *fqplane[17];
41    //TMatrixD *fnqplane[17];
42    //TMatrixF *fmatrix[17];
43    //TMatrixF *fnmat[17];
44    //TMatrixD *fmatrix;
45    //TMatrixD *fnmat;
46    TMatrixF *fmatrix;
47    //TMatrixF *fnmat;
48    TMatrixF *fnmat[17];
49    //Int_t finmat[43][191];
50    
51  //===============================================================================  //===============================================================================
52  bool Select( PamLevel2* event ){  bool Select( PamLevel2* event ){
# Line 40  bool Select( PamLevel2* event ){ Line 62  bool Select( PamLevel2* event ){
62    // tracker pre-selection    // tracker pre-selection
63    //------------------------------------------------------------------    //------------------------------------------------------------------
64    TrkTrack *trk = track->GetTrkTrack();    TrkTrack *trk = track->GetTrkTrack();
65      float rigidity = trk->GetRigidity();
66      if ( CRIG ) rigidity = event->GetCaloLevel2()->qtot/260.;
67      if ( SRIG ) rigidity = event->GetGPamela()->P0;
68    bool TRACK__OK = false;    bool TRACK__OK = false;
69    if(    if(
70       trk->chi2 >0     &&       trk->chi2 >0     &&
# Line 66  bool Select( PamLevel2* event ){ Line 90  bool Select( PamLevel2* event ){
90       !event->GetAcLevel2()->CARDhit() &&       !event->GetAcLevel2()->CARDhit() &&
91       !event->GetAcLevel2()->CAThit() &&       !event->GetAcLevel2()->CAThit() &&
92       true ) TOF__OK = true;       true ) TOF__OK = true;
93    if( !TOF__OK )return false;    if( !TOF__OK && !SIMU)return false;
94    //------------------------------------------------------    //------------------------------------------------------
95    // no albedo    // no albedo
96    //------------------------------------------------------    //------------------------------------------------------
97    if(    track->GetToFTrack()->beta[12]<=0.2 ||    if(    !SIMU && (track->GetToFTrack()->beta[12]<=0.2 ||
98           track->GetToFTrack()->beta[12] >= 1.5 ) return false;                     track->GetToFTrack()->beta[12] >= 1.5) ) return false;
99    
100    //------------------------------------------------------    //------------------------------------------------------
101    bool CUT1 = false;    bool CUT1 = false;
102    if(    if(
103       trk->nstep<100   &&       trk->nstep<100   &&
104       trk->GetRigidity()<400.   &&       rigidity<400.   &&
105       trk->GetRigidity()>0.1   &&       rigidity>0.1   &&
      trk->GetDeflection()<0.   &&  
106       trk->resx[0]<0.001 &&       trk->resx[0]<0.001 &&
107       trk->resx[5]<0.001 &&       trk->resx[5]<0.001 &&
108       track->IsSolved() &&       track->IsSolved() &&
# Line 87  bool Select( PamLevel2* event ){ Line 110  bool Select( PamLevel2* event ){
110       true ) CUT1 = true;       true ) CUT1 = true;
111    //------------------------------------------------------    //------------------------------------------------------
112    if( !CUT1 )return false;    if( !CUT1 )return false;
113      if ( trk->GetDeflection()>0. && !SIMU ) return false;
114    
115      //
116      //  ELENA'S CUT
117      //
118      //
119      // lever-arm 6
120      //====================================================
121      bool LX6=false;
122      if(
123         track->GetTrkTrack()->GetLeverArmX()==6 &&
124         !track->GetTrkTrack()->IsBad(0,0) &&
125         !track->GetTrkTrack()->IsBad(5,0) &&
126         track->GetTrkTrack()->resx[0]<0.001 &&
127         track->GetTrkTrack()->resx[5]<0.001 &&
128         track->GetTrkTrack()->IsInsideCavity() &&
129         true ) LX6 = true;
130      
131      //====================================================
132      // lever-arm 5
133      //====================================================
134      bool LX5=false;
135      if(
136         track->GetTrkTrack()->GetLeverArmX()==5 &&
137         true ){
138        if(
139           track->GetTrkTrack()->XGood(0) && track->GetTrkTrack()->XGood(4)
140           ){
141          
142          if(
143             !track->GetTrkTrack()->IsBad(0,0) &&
144             !track->GetTrkTrack()->IsBad(4,0) &&
145             track->GetTrkTrack()->resx[0]<0.001 &&
146             track->GetTrkTrack()->resx[4]<0.001 &&
147             track->GetTrkTrack()->IsInsideCavity() &&
148             true) LX5 = true;        
149        }else if (
150                  track->GetTrkTrack()->XGood(1) && track->GetTrkTrack()->XGood(5)
151                  ){
152          
153          if(
154             !track->GetTrkTrack()->IsBad(1,0) &&
155             !track->GetTrkTrack()->IsBad(5,0) &&
156             track->GetTrkTrack()->resx[1]<0.001 &&
157             track->GetTrkTrack()->resx[5]<0.001 &&
158             track->GetTrkTrack()->IsInsideCavity() &&
159             true) LX5 = true;
160        }
161      }
162      if ( !LX5 && !LX6 ) return false;
163    Float_t defl = trk->GetDeflection();    Float_t defl = trk->GetDeflection();
164    float p0 = 1.111588e+00;    float p0 = 1.111588e+00;
165    float p1 = 1.707656e+00;    float p1 = 1.707656e+00;
# Line 104  bool Select( PamLevel2* event ){ Line 177  bool Select( PamLevel2* event ){
177       track->GetTrkTrack()->chi2 < chi2m &&       track->GetTrkTrack()->chi2 < chi2m &&
178       true ) CUT2 = true;         true ) CUT2 = true;  
179    if ( !CUT2 ) return false;    if ( !CUT2 ) return false;
180      float dedxtrk    = trk->GetDEDX();
181      //  float zcutn =  9.   + 20./(rigidity*rigidity);
182      float zcut2 =  3.   + 4.3/(rigidity*rigidity);
183      float zcut1 =  0.52 + 0.455/(rigidity*rigidity);
184      Bool_t Z1 = false;
185      if(dedxtrk > zcut1 && dedxtrk < zcut2){
186        Z1=true;
187      }
188      if ( !Z1 && !SIMU ) return false;
189    //------------------------------------------------------    //------------------------------------------------------
190    //    //
191    // energy momentum match    // energy momentum match
# Line 115  bool Select( PamLevel2* event ){ Line 197  bool Select( PamLevel2* event ){
197    //    //
198    for (Int_t i=0; i < 22; i++){    for (Int_t i=0; i < 22; i++){
199      if ( track->GetCaloTrack()->tibar[i][1] < 0 || track->GetCaloTrack()->tibar[i][0] < 0 ){      if ( track->GetCaloTrack()->tibar[i][1] < 0 || track->GetCaloTrack()->tibar[i][0] < 0 ){
200          return false;        return false;
201        };
202      };
203      //
204      if ( event->GetCaloLevel2()->qtot == 0. ) return false;
205      if ( rigidity>5. && track->GetCaloTrack()->qtrack/event->GetCaloLevel2()->qtot < 0.4 ) return false;
206      if ( rigidity<1. && track->GetToFTrack()->beta[12] < 0.8 ) return false;
207      if ( rigidity>50. ){
208        if ( trk->GetNX()<5  &&
209             trk->GetNY()<4 ) return false;
210        //
211        Bool_t sphit = false;
212        for ( Int_t plane = 0; plane < 6; plane++){
213          if ( !trk->XGood(plane) ){
214            for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsx(); sing++){
215              TClonesArray &t = *(event->GetTrkLevel2()->SingletX);
216              TrkSinglet *singlet = (TrkSinglet*)t[sing];
217              if ( (singlet->plane-1) == plane ){
218                Float_t x = (singlet->coord[0]+singlet->coord[1])/2.;
219                if ( fabs(track->GetTrkTrack()->xv[plane] - x) < 1. ) sphit = true;
220              };
221            };
222          };
223          if ( !trk->YGood(plane) ){
224            for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsy(); sing++){
225              TClonesArray &t = *(event->GetTrkLevel2()->SingletY);
226              TrkSinglet *singlet = (TrkSinglet*)t[sing];
227              if ( (singlet->plane-1) == plane ){
228                Float_t x1 = (singlet->coord[0]);
229                Float_t x2 = (singlet->coord[1]);
230                if ( fabs(track->GetTrkTrack()->yv[plane] - x1) < 1. ) sphit = true;
231                if ( fabs(track->GetTrkTrack()->yv[plane] - x2) < 1. ) sphit = true;
232              };
233            };
234          };
235      };      };
236        if ( sphit ) return false; // spurious hit along the track
237    };    };
238    //    //
239    if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;    Int_t ti0 = track->GetCaloTrack()->tibar[0][1]-1;
240      
241      Int_t view = 0;
242      Int_t plane = 0;
243      Int_t strip = 0;
244      Float_t mip = 0.;
245      //
246      for ( Int_t i=0; i<event->GetCaloLevel1()->istrip; i++ ){
247        //
248        mip = event->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
249        if ( view == 1 && plane == 0 && strip == ti0 && mip > 4.) return false;
250        if ( view == 1 && (plane >0 || strip > ti0) ) break;
251      };
252      //  if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;
253    
254      if ( rigidity > 2.2 || rigidity < 1.5 ) return false;
255      //  printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG);
256    //    //
257    return true;    return true;
258  }  }
# Line 134  bool Select( PamLevel2* event ){ Line 267  bool Select( PamLevel2* event ){
267  void CreateHistos( PamLevel2* event , TString file){  void CreateHistos( PamLevel2* event , TString file){
268    
269    cf = new CaloFranzini(event);    cf = new CaloFranzini(event);
270    cf->CreateMatrixFile(file.Data());    //
271      if ( MATRIX ){
272        cf->UpdateMatrixFile(file.Data());
273        cf->LoadBin();
274        if ( !FULL ){
275          cf->LoadLong();
276        } else {
277          cf->LoadFull();
278        };
279      } else {
280        cf->CreateMatrixFile(file.Data());
281      };
282      //
283    //    //
284    nbin = 18;    nbin = 18;
285    rig[0] = 0.1;    rig[0] = 0.1;
# Line 156  void CreateHistos( PamLevel2* event , TS Line 301  void CreateHistos( PamLevel2* event , TS
301    rig[16] = 200.;    rig[16] = 200.;
302    rig[17] = 400.;    rig[17] = 400.;
303    //    //
304  //  memset(qplane, 0, 44*18*sizeof(Float_t));    memset(rmean, 0, 17*sizeof(Float_t));
305    memset(rmean, 0, 18*sizeof(Float_t));    memset(ntot, 0, 17*sizeof(Int_t));
306    memset(qqplane, 0, 43*18*sizeof(Float_t));    //  memset(finmat, 0, 43*191*sizeof(Int_t));  
307    memset(nnqplane, 0, 43*18*sizeof(Float_t));    //
308    //  //  for (Int_t i=0; i < 17 ; i++){
309    for (Int_t i=0; i < 18 ; i++){      for (Int_t i=3; i < 4 ; i++){
310      matrix[i] = new TMatrixD(43,43);      if ( !FULL ){
311      qplane[i] = new TArrayF(43);        matrix[i] = new TMatrixD(43,43);
312      nqplane[i] = new TArrayI(43);        qplane[i] = new TArrayF(43);
313      nmat[i] = new TMatrixD(43,43);        nqplane[i] = new TArrayI(43);
314    };        nmat[i] = new TMatrixD(43,43);
315  //  for (Int_t i=0; i < 18 ; i++){      } else {
316  //      for (Int_t j = 0; j < 43;i++){        if ( MATRIX ){
317  //              (*qplane[i])[j] = 0.;          //      fmatrix = new TMatrixF(4128,4128);
318  //              qplane[i]->Reset();          //      fnmat = new TMatrixF(4128,4128);
319  //              nqplane[i]->Reset();          //      fmatrix = new TMatrixF(8213,8213);
320  //              for (Int_t k = 0; k < 43;k++){          //      fnmat = new TMatrixF(8213,8213);
321  //                      (*matrix[i])[j][k] = 0.;          fmatrix = new TMatrixF(MDIM,MDIM);
322  //                      (*nmat[i])[j][k] = 0.;          //      fnmat = new TMatrixF(MDIM,MDIM);
323  //              };          fnmat[i] = new TMatrixF(43,191);
324  //      };  //      cf->WriteFullMatrix(fmatrix, i);
325  //  };          //      cf->WriteFullNMatrix(fnmat, i);
326    //      delete fmatrix;
327            //      delete fnmat;
328            //fnmat[i] = new TMatrixI(8213,8213);
329          } else {
330            fqplane = new TMatrixD(43,191); // 43 planes x 191 strip (= 1 + 95 x 2, one strip is the one transversed by the track that could be on the extreme right or left)
331            fnqplane = new TMatrixD(43,191);//
332            //
333            cf->WriteFullMean(fqplane, i);
334            cf->WriteFullNMean(fnqplane, i);
335            delete fqplane;
336            delete fnqplane;
337            //
338          };
339        };
340      };
341    //    //
342  }  }
343    
344  //===============================================================  //===============================================================
345  void FindAverage( PamLevel2* L2, int iev ){  void FindAverage( PamLevel2* L2, int iev ){
346    //    //
   //  printf("SELEZIONATO \n");  
347    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
348      if ( SRIG ) erig = L2->GetGPamela()->P0;
349      if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.;
350    //    //
351    Int_t rbi = 0;    Int_t rbi = 0;
352    for (Int_t i = 0; i<nbin-1; i++){ // < 17    for (Int_t i = 0; i<nbin-1; i++){
353      if ( erig>=rig[i] && erig < rig[i+1] ){      if ( erig>=rig[i] && erig < rig[i+1] ){
354        rbi = i;        rbi = i;
355        break;        break;
356      };      };
357    };    };
358    //    //
   //  printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]);  
   //  
359    if ( erig < rig[0] ) return;    if ( erig < rig[0] ) return;
360    if ( erig >= rig[nbin-1] ) return;    if ( erig >= rig[nbin-1] ) return;
361    //    //
362    rmean[rbi] += erig;    rmean[rbi] += erig;
363      ntot[rbi]++;
364    //    //
365    Int_t dgf = 43;    if (!FULL ){
366    //      Int_t dgf = 43;
367  //  for (Int_t i=0; i < 22; i++){      //
368  //    if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){      for (Int_t i=0; i<dgf; i++){
369  //      dgf = 2 * i;        (*nqplane[rbi])[i]++;
370   //     break;      };
371    //  };      //
372   //   if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){      // Fill the estrip matrix
373   //     dgf = 1 + 2 * i;      //
374   //     break;      Int_t nplane = 0;
375   //   };      Int_t view = 0;
376   // };      Int_t plane = 0;
377    //      Int_t strip = 0;
378  //  if ( dgf < 43 && dgf > 37 ) dgf--;      Float_t mip = 0.;
379    //      //
380    for (Int_t i=0; i<dgf; i++){      for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
381      (*nqplane[rbi])[i]++;        //
382          nnqplane[rbi][i]++;        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
383    };        //
384    //        nplane = 1 - view + 2 * plane;
385    //  printf(" SELEZIONATO rig %f dgf %i \n",erig,dgf);        if ( erig > 4. && nplane == 0 && mip > 15. ) printf(" IEV %i erig %f OBT %u pkt %u file %s \n",iev,erig,L2->GetOrbitalInfo()->OBT,L2->GetOrbitalInfo()->pkt_num,L2->GetPamTree()->GetFile()->GetName());
386    //        //printf(" IEV %i OBT %u pkt %u file %s \n",iev,L2->GetOrbitalInfo()->OBT,L2->GetOrbitalInfo()->pkt_num,L2->GetPamTree()->GetFile()->GetName());
387    //        if ( nplane > 37 ) nplane--;
388    // Fill the estrip matrix        if ( nplane < dgf ){
389    //          (*qplane[rbi])[nplane] += mip;
390    Int_t nplane = 0;        };
391    Int_t view = 0;        //
392    Int_t plane = 0;      };
393    Int_t strip = 0;    } else {
394    Float_t mip = 0.;      //
395    //      // FULL CALORIMETER
   for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){  
396      //      //
397      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      fqplane = cf->LoadFullAverage(rbi);
398        fnqplane = cf->LoadFullNAverage(rbi);
399        CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
400      //      //
401      nplane = 1 - view + 2 * plane;      Int_t nplane = 0;
402      if ( nplane > 37 ) nplane--;      Int_t view = 0;
403      //    printf(" mip %f view %i plane %i => nplane %i \n",mip,view,plane,nplane);      Int_t plane = 0;
404      if ( nplane < dgf ){      Int_t strip = 0;
405        //      printf(" PRIMA: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]);      Float_t mip = 0.;
406        (*qplane[rbi])[nplane] += mip;      //
407          qqplane[rbi][nplane] += mip;      Int_t cs = 0;
408        //      printf(" DOPO: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]);      Int_t cd = 0;
409        Int_t mstrip = 0;
410        //
411        for (Int_t j=0; j<2; j++){
412          for (Int_t i=0; i<21; i++){
413            nplane = 1 - j + 2*i;
414            if ( nplane > 37 ) nplane--;
415            //
416            cs = ct->tibar[i][j] - 1;
417            //
418            cd = 95 - cs;  
419            //
420            for (Int_t k=0; k<191; k++){
421              mstrip = cd + k;
422              //      if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.;
423              if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.;
424            };
425          };
426      };      };
427      //      //
428        //
429        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
430          //
431          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
432          //
433          nplane = 1 - view + 2 * plane;
434          if ( nplane > 37 ) nplane--;
435          //
436          cs = ct->tibar[plane][view] - 1;
437          //
438          cd = 95 - cs;
439          //
440          mstrip = cd + strip;
441          //
442          //      (*fqplane[rbi])[nplane][mstrip] += mip;
443          (*fqplane)[nplane][mstrip] += mip;
444          //
445        };
446        //
447        cf->WriteFullMean(fqplane, rbi);
448        cf->WriteFullNMean(fnqplane, rbi);
449        cf->UnLoadFullAverage(rbi);
450        cf->UnLoadFullNAverage(rbi);
451        delete fqplane;
452        delete fnqplane;
453        //
454    };    };
455  }  }
456    
457    void CalculateAverage(){
458      //  
459      if ( !FULL ){
460        for (Int_t i=0; i<nbin-1; i++){
461          if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];
462          for (Int_t j=0; j<43 ; j++){
463            if ( (*nqplane[i])[j] > 0 ){  
464              (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];
465            } else {
466              (*qplane[i])[j] = 0.;
467            };
468            printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]);
469          };
470        };
471        for (Int_t i=0; i<nbin-1; i++){
472          //
473          cf->WriteLongMean(qplane[i], i);
474          //
475        };
476      } else {
477        //
478        for (Int_t i=0; i<nbin-1; i++){
479          fqplane = cf->LoadFullAverage(i);
480          fnqplane = cf->LoadFullNAverage(i);
481          if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]);
482          //
483          for (Int_t j=0; j<43 ; j++){
484            for (Int_t k=0; k<191; k++){
485    //        if ( (*fnqplane[i])[j][k] > 0 ){  
486    //          (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
487    //        } else {
488    //          (*fqplane[i])[j][k] = 0.;
489    //        };
490    //        printf(" BIN %i plane %i strip %i average energy %f qplane %f nqplane %f \n",i,j,k,rmean[i],(*fqplane[i])[j][k],(*fnqplane[i])[j][k]);
491              if ( (*fnqplane)[j][k] > 0 ){  
492                if ( (*fqplane)[j][k] == 0. ) (*fqplane)[j][k] = 0.7;
493                (*fqplane)[j][k] /= (Float_t)(*fnqplane)[j][k];
494              } else {
495                (*fqplane)[j][k] = 0.;
496              };
497              printf(" BIN %i plane %i strip %i average energy %f qplane %f nqplane %f \n",i,j,k,rmean[i],(*fqplane)[j][k],(*fnqplane)[j][k]);
498            };
499          };
500          cf->WriteFullMean(fqplane, i);
501          cf->WriteFullNMean(fnqplane, i);
502          cf->UnLoadFullAverage(i);
503          cf->UnLoadFullNAverage(i);
504          delete fqplane;
505          delete fnqplane;
506        };
507        //
508    //     for (Int_t i=0; i<nbin-1; i++){
509    //       //
510    //       cf->WriteFullMean(fqplane[i], i);
511    //       //
512    //     };
513      };
514      //  
515      cf->WriteNumBin(nbin);
516      //
517      TArrayF *rigbin = new TArrayF(18, rig);
518      cf->WriteRigBin(rigbin);
519      //
520      TArrayF *rmeanbin = new TArrayF(17, rmean);
521      TFile *file =  cf->GetFile();
522      file->cd();
523      file->WriteObject(&(*rmeanbin), "binrigmean");
524      //
525      //
526    }
527    
528  //===============================================================  //===============================================================
529  void FindMatrix( PamLevel2* L2, int iev ){  void FindMatrix( PamLevel2* L2, int iev ){
530    //    //
   //  printf("SELEZIONATO \n");  
531    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
532      if ( SRIG ) erig = L2->GetGPamela()->P0;
533      if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.;
534    //    //
535    Int_t rbi = 0;    Int_t rbi = 0;
536    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
# Line 264  void FindMatrix( PamLevel2* L2, int iev Line 540  void FindMatrix( PamLevel2* L2, int iev
540      };      };
541    };    };
542    //    //
543    //  printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]);  if ( rbi != 3 ) return;
   //  
544    if ( erig < rig[0] ) return;    if ( erig < rig[0] ) return;
545    if ( erig >= rig[nbin-1] ) return;    if ( erig >= rig[nbin-1] ) return;
546    //    //
547    Int_t dgf = 43;    if ( !FULL ){
548    //      Int_t dgf = 43;
549  //  for (Int_t i=0; i < 22; i++){      //
550  //    if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){      for (Int_t i=0; i<dgf; i++){
551  //      dgf = 2 * i;        for (Int_t j=0; j<dgf; j++){
552   //     break;          (*nmat[rbi])[i][j] += 1.;
553   //   };        };
  //   if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){  
  //     dgf = 1 + 2 * i;  
  //     break;  
  //   };  
 //  };  
   //  
 //  if ( dgf < 43 && dgf > 37 ) dgf--;  
   //  
   for (Int_t i=0; i<dgf; i++){  
     for (Int_t j=0; j<dgf; j++){  
       (*nmat[rbi])[i][j] += 1.;  
554      };      };
   };  
   //  
   //  printf(" SELEZIONATO rig %f dgf %i \n",erig,dgf);  
   //  
   //  
   // Fill the estrip matrix  
   //  
   Int_t nplane = 0;  
   Int_t view = 0;  
   Int_t plane = 0;  
   Int_t strip = 0;  
   Float_t mip = 0.;  
   Float_t hpl[43];  
   memset(hpl,0,43*sizeof(Float_t));  
   for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){  
555      //      //
556      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      // Fill the estrip matrix
557      //      //
558      nplane = 1 - view + 2 * plane;      Int_t nplane = 0;
559      if ( nplane > 37 ) nplane--;      Int_t view = 0;
560      if ( nplane < dgf ){      Int_t plane = 0;
561        hpl[nplane] += mip;      Int_t strip = 0;
562        Float_t mip = 0.;
563        Float_t hpl[43];
564        memset(hpl,0,43*sizeof(Float_t));
565        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
566          //
567          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
568          //
569          nplane = 1 - view + 2 * plane;
570          if ( nplane > 37 ) nplane--;
571          if ( nplane < dgf ){
572            hpl[nplane] += mip;
573          };
574          //
575      };      };
576      //      //
577    };      for (Int_t i=0; i<dgf; i++){
578    //        for (Int_t j=0; j<dgf; j++){
579    for (Int_t i=0; i<dgf; i++){          (*matrix[rbi])[i][j] += (hpl[i] - cf->GetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig));
580      for (Int_t j=0; j<dgf; j++){        };
       // if ( rbi == 1 ) printf(" PRIMA: matrix %f qpli %f qplj %f hpli %f hplj %f \n",(*matrix[rbi])[i][j],(*qplane[rbi])[i],(*qplane[rbi])[j],hpl[i],hpl[j]);  
       (*matrix[rbi])[i][j] += (hpl[i] - (*qplane[rbi])[i]) * (hpl[j] - (*qplane[rbi])[j]);  
       //      if ( rbi == 1 ) printf(" DOPO: matrix %f qpli %f qplj %f hpli %f hplj %f \n",(*matrix[rbi])[i][j],(*qplane[rbi])[i],(*qplane[rbi])[j],hpl[i],hpl[j]);  
581      };      };
582    };    } else {
583  }      //
584        // FULL CALORIMETER
585  void CalculateAverage(){      //
586    for (Int_t i=0; i<nbin; i++){      printf(" Retrieve matrix %i IEV %i \n",rbi,iev);
587      if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];      //    fmatrix = cf->LoadFullMatrix(rbi);
588      for (Int_t j=0; j<43 ; j++){  //    cf->LoadFullMatrix(rbi,fmatrix);
589        if ( (*nqplane[i])[j] > 0 ){        //    fnmat = cf->LoadFullNMatrix(rbi);
590          (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];      printf(" done \n");
591        } else {      printf(" start loop \n");
592          (*qplane[i])[j] = 0.;      //
593        CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
594        //
595        Int_t nplane = 0;
596        Int_t view = 0;
597        Int_t plane = 0;
598        Int_t strip = 0;
599        Float_t mip = 0.;
600        //
601        //    Int_t mindgf = 48;
602        //    Int_t dgf = 143;
603    //       Int_t mindgf = 0; //tutto
604    //        Int_t dgf = 191; //tutto
605    //    Int_t mindgf = 94;
606    //    Int_t dgf = 96;
607        Int_t mindgf = 84;
608        Int_t dgf = 106;
609        Int_t cs = 0;
610        Int_t cd = 0;
611        Int_t mstrip = 0;
612        //
613        Float_t mipv[43][191];
614        memset(mipv,0,43*191*sizeof(Float_t));
615        //
616        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
617          //
618          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
619          //
620          nplane = 1 - view + 2 * plane;
621          if ( nplane > 37 ) nplane--;
622          //
623          cs = ct->tibar[plane][view] - 1;
624          //
625          cd = 95 - cs;
626          //
627          mstrip = cd + strip;
628          //
629          mipv[nplane][mstrip] = mip;
630          //
631        };
632        //
633        Float_t mip1 = 1.;
634        Int_t cs1;
635        Int_t cd1;
636        Float_t mip2 = 1.;
637        Int_t cs2;
638        Int_t cd2;
639        Int_t mi = -1;
640        Int_t mj = -1;
641        Int_t nn1 = 0;
642        Int_t pl1 = 0;
643        Int_t vi1 = 0;
644        Int_t nn2 = 0;
645        Int_t pl2 = 0;
646        Int_t vi2 = 0;
647        Int_t mstrip1min = 0;
648        Int_t mstrip1max = 0;
649        Int_t mstrip2min = 0;
650        Int_t mstrip2max = 0;
651        //
652        Int_t toto = 0;
653        //
654        for (Int_t nplane1=0; nplane1<43; nplane1++){
655          if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
656          vi1 = 1;
657          if ( nn1%2 ) vi1 = 0;
658          pl1 = (nn1 - 1 + vi1)/2;
659          //
660          cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1
661          //
662          cd1 = 95 - cs1;
663          //
664          mstrip1min = TMath::Max(mindgf,(cd1+0));
665          mstrip1max = TMath::Min(dgf,(cd1+95)) + 1;
666          //
667          if ( nplane1 == 0 || nplane1 == 42 ) printf(" pl %i mstrip1min %i mstrip1max %i \n",nplane1,mstrip1min,mstrip1max);
668          //
669          for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
670            //      printf(".\n");
671            //
672            mj = -1;
673            //
674            mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);
675            //
676            mi = (nplane1 * 191) + mstrip1;
677            //
678            //      if ( mstrip1 > mstrip1min ) break;
679            //      if ( mstrip1 > dgf ) break;
680            //      if ( mstrip1 >= mindgf && mstrip1 <= dgf && mstrip1 >= mstrip1min && mstrip1 <= mstrip1max ){
681            //
682            //      finmat[nplane1][mstrip1]++;
683            (*fnmat[rbi])[nplane1][mstrip1] += 1.;
684            //
685            if ( mip1 != 0. ){
686              //
687              for (Int_t nplane2=0; nplane2<43; nplane2++){
688                //
689                if ( nplane2 >= 37 ) nn2 = nplane2 + 1;
690                vi2 = 1;
691                if ( nn2%2 ) vi2 = 0;
692                pl1 = (nn2 - 1 + vi2)/2;
693                //
694                cs2 = ct->tibar[pl2][vi2] - 1;
695                //
696                cd2 = 95 - cs2;
697                //
698                //    mstrip2min = cd2 + 0;
699                //    mstrip2max = cd2 + 95;
700                mstrip2min = TMath::Max(mindgf,(cd2+0));
701                mstrip2max = TMath::Min(dgf,(cd2+95)) + 1;
702                //
703                if ( nplane1 == 0 && nplane2 == 0 && mstrip1==mstrip1min ) printf(" mstrip2min %i mstrip2max %i \n",mstrip2min,mstrip2max);
704                //
705                for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
706                  //
707                  mip2 = mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi);
708                  //
709                  if ( mip2 != 0. ){
710                    //
711                    mj = (nplane2 * 191) + mstrip2;
712                    //            mj++;
713                    //
714                    //            if ( mstrip2 > mstrip2min ) break;
715                    //            if ( mstrip2 > dgf ) break;
716                    //            if ( mstrip2 >= mindgf && mstrip2 <= dgf && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max ){
717                    //              if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){
718                    //              (*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig));
719                    //                (*fnmat[rbi])[mi][mj] += 1.;
720                    (*fmatrix)[mi][mj] += (mip1 * mip2); // giusto
721    //              (*fmatrix)[mi][mj] += (mip1 * mip2) * 1000000.;
722                    toto++;
723                    //            (*fmatrix)[mi][mj] += 1.;
724                    //              cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);
725                    //              cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi);
726                    //                (*fnmat)[mi][mj] += 1.;
727                    //              };
728                  };
729                };
730              };
731            };
732        };        };
       printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]);  
 //      if ( nnqplane[i][j] > 0 ){    
 //      qqplane[i][j] /= (Float_t)nnqplane[i][j];  
 //      } else {  
 //      qqplane[i][j] = 0.;  
 //      };  
 //      printf(" BBIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],qqplane[i][j],nnqplane[i][j]);  
733      };      };
734    };      //  
735          printf(" toto = %i \n",toto);
736    cf->WriteNumBin(nbin);      printf("\n done \n");
737        printf(" write matrix \n");
738    TArrayF *rigbin = new TArrayF(18, rig);  //    cf->WriteFullMatrix(fmatrix, rbi);
739    cf->WriteRigBin(rigbin);      //    cf->WriteFullNMatrix(fnmat, rbi);
740        printf(" done \n");
741    TArrayF *rmeanbin = new TArrayF(18, rmean);      printf(" unload matrix \n");
742    TFile *file =  cf->GetFile();   //   cf->UnLoadFullMatrix(rbi);
743    file->cd();      //    cf->UnLoadFullNMatrix(rbi);
744    //  rigbin->Write("binrig");      printf(" done \n");
745    file->WriteObject(&(*rmeanbin), "binrigmean");      printf(" delete matrix \n");
746    //    delete fmatrix;
747  //  cf->WriteRigBin(rigbin);      //    delete fnmat;
748          printf(" done \n");
   for (Int_t i=0; i<nbin; i++){  
       
     cf->WriteLongMean(qplane[i], i);  
   
749    };    };
750  }  }
751    
752  //===============================================================  //===============================================================
753  // Save histograms  // Save histograms
754  //  //
# Line 370  void CalculateAverage(){ Line 758  void CalculateAverage(){
758  //  //
759  //===============================================================  //===============================================================
760  void SaveHistos(){  void SaveHistos(){
761        //
762    printf("Finished, calculating average and inverting matrices\n");    if ( MATRIX ){
   
   
     
 //  cf->WriteNumBin(nbin);  
   
 //  TArrayF *rigbin = new TArrayF(18, rig);  
 //  cf->WriteRigBin(rigbin);  
     
   for (Int_t i=0; i<nbin; i++){  
       
 //    cf->WriteLongMean(qplane[i], i);  
   
763      //      //
764      // determine the average matrix      printf("Finished, calculating average and inverting matrices\n");
765      //          //
766      for (Int_t ii=0; ii<43; ii++){      if ( !FULL ){
767        for (Int_t j=0; j<43; j++){        for (Int_t i=0; i<nbin-1; i++){
768          //      if ( i == 1 ) printf(" num %f \n",(*nmat[i])[ii][j]);          //
769          if ( (*nmat[i])[ii][j] > 0. ){          // determine the average matrix
770            //      if ( i == 1 ) printf(" value %f \n",(*matrix[i])[ii][j]);          //    
771            (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];          for (Int_t ii=0; ii<43; ii++){
772              for (Int_t j=0; j<43; j++){
773                if ( (*nmat[i])[ii][j] > 0. ){
774                  (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];
775                } else {
776                  (*matrix[i])[ii][j] = 0.;
777                };
778              };
779            };
780            //
781            cf->WriteLongMatrix(matrix[i],i);
782            //
783            if ( matrix[i]->Determinant() == 0. ){
784              printf("\n");
785              for (Int_t ii=0; ii<43; ii++){
786                for (Int_t j=0; j<43; j++){
787                  printf(" %.f",(*matrix[i])[ii][j]);        
788                };
789                printf("\n");
790              };
791              printf("\n");
792              printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
793          } else {          } else {
794            (*matrix[i])[ii][j] = 0.;            Double_t det = 0.;
795              TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));
796              printf(" Bin %i determinant is %f \n",i,det);
797              cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);
798          };          };
799        };        };
800      };      } else {
801              //
802  //     if ( i == 2 ){        // FULL
803  //       printf("\n");        //
804  //       for (Int_t ii=0; ii<43; ii++){        //      for (Int_t i=0; i<nbin-1; i++){
805  //      for (Int_t j=0; j<43; j++){  //      for (Int_t i=3; i<5; i++){
806  //        printf(" %.f",(*matrix[i])[ii][j]);            for (Int_t i=3; i<4; i++){
807            //
808            // determine the average matrix
809            //    
810    //      fmatrix = cf->LoadFullMatrix(i);
811            //      fnmat = cf->LoadFullNMatrix(i);
812            //
813    //      for (Int_t ii=0; ii<MDIM; ii++){
814    //        for (Int_t j=0; j<MDIM; j++){
815    // //       if ( (*fnmat[i])[ii][j] > 0. ){
816    // //         (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
817    // //       } else {
818    // //         (*fmatrix[i])[ii][j] = 0.;
819    // //       };
820    //          if ( (*fnmat)[ii][j] > 0. ){
821    //            (*fmatrix)[ii][j] /= (*fnmat)[ii][j];
822    //          } else {
823    //            (*fmatrix)[ii][j] = 0.;
824    //          };
825    //        };
826  //      };  //      };
827  //      printf("\n");  //
828  //       };  //      TMatrixF *mymat = new TMatrixF(129,129);
829  //       printf("\n");          TMatrixF *mymat = new TMatrixF(989,989);
830  //     };          Int_t i1 = -1;
831                Int_t j1 = -1;
832            int mi,mj;
833            Int_t nonzero = 0;
834            Int_t nonzero1 = 0;
835            for (Int_t ii=0; ii<43; ii++){
836              for (Int_t j=0; j<191; j++){
837    //          if ( (*fnmat[i])[ii][j] > 0. ){
838    //            (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
839    //          } else {
840    //            (*fmatrix[i])[ii][j] = 0.;
841    //          };
842                i1 =  (ii * 191) + j;
843                //      j1 = -1;
844                for (Int_t iij=0; iij<43; iij++){
845                  for (Int_t jj=0; jj<191; jj++){
846                    //
847                    j1 = (iij * 191) + jj;
848                    //              j1++;
849                    //              if ( finmat[ii][j] > 0 ){
850                    //                (*fmatrix)[i1][j1] /= finmat[ii][j];
851                    if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1] == (*fmatrix)[i1][j1]) ){
852                      (*fmatrix)[i1][j1] = 1.;
853                    } else {
854                      (*fmatrix)[i1][j1] /= (*fnmat[i])[ii][j];
855                      nonzero++;
856                      if ( i1 == 0 ) nonzero1++;
857                    };
858                    //              if ( j>=94 && j <=96 && jj >=94 && jj<=96 ){
859                    //                      mi = (ii*3) + j -94;
860                    //                      mj = (iij*3) + jj -94;
861                    //                      (*mymat)[mi][mj] = (*fmatrix)[i1][j1];
862                    //              };
863    
864    
865                    if ( j>=84 && j <=106 && jj >=84 && jj<=106 ){
866                            mi = (ii*3) + j -84;
867                            mj = (iij*3) + jj -84;
868                            (*mymat)[mi][mj] = (*fmatrix)[i1][j1];
869                    };
870    
871                  };
872                };
873              };
874            };
875            //
876            printf(" Matrix has %i non-zero elements \n",nonzero);  
877            printf(" Matrix has %i non-zero elements on the first row\n",nonzero1);
878            //
879    //      Bool_t BAD = false;
880    //      for (Int_t ii=0; ii<43; ii++){
881    //        for (Int_t j=0; j<191; j++){
882    //          //
883    //          i1 =  (ii * 191) + j;
884    //          //
885    //          for (Int_t iij=0; iij<43; iij++){
886    //            for (Int_t jj=0; jj<191; jj++){
887    //              //
888    //              j1 = (iij * 191) + jj;
889    //              //
890    //              //              printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
891    //              if ( (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1]==(*fmatrix)[i1][j1]) ){
892    //                printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
893    //                printf(" che schifo! \n");
894    //                BAD = true;
895    //              };
896    //              //
897    //            };
898    //          };
899    //        };
900    //      };
901    //      //
902    //      if ( BAD ) printf(" questa matrice fa cagare \n");
903            //
904            //
905            //      cf->WriteFullMatrix(fmatrix[i],i);
906            cf->WriteFullMatrix(fmatrix, i);
907            //      cf->WriteFullNMatrix(fnmat, i);
908            cf->WriteFullNMatrix(fnmat[i], i);
909            //
910            //      if ( fmatrix[i]->Determinant() == 0. ){
911            if ( fmatrix->Determinant() == 0. ){
912              printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
913            } else {
914              Double_t det = 0.;
915              //      TMatrixF invmatrix = (TMatrixF)(fmatrix[i]->Invert(&det));
916              //      printf(" Bin %i determinant is %f \n",i,det);
917              //      cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);
918              TMatrixF invmatrix = (TMatrixF)(fmatrix->Invert(&det));
919              printf(" Bin %i determinant is %f \n",i,det);
920              cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);
921            };
922    
923            if ( mymat->Determinant() == 0. ){
924              printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
925            } else {
926              Double_t det = 0.;
927              TMatrixF invmatrix = (TMatrixF)(mymat->Invert(&det));
928              printf(" Bin %i determinant is %f \n",i,det);
929              cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);
930            };
931            cf->WriteFullMatrix(mymat, 99);
932    
933      cf->WriteLongMatrix(matrix[i],i);  
934                //
935      if ( matrix[i]->Determinant() == 0. ){          cf->UnLoadFullMatrix(i);
936        printf("\n");          //      cf->UnLoadFullNMatrix(i);
937        for (Int_t ii=0; ii<43; ii++){          delete fmatrix;
938          for (Int_t j=0; j<43; j++){          //      delete fnmat;
939            printf(" %.f",(*matrix[i])[ii][j]);              //
         };  
         printf("\n");  
940        };        };
       printf("\n");  
       printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);  
     } else {  
       Double_t det = 0.;  
       TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));  
       printf(" Bin %i determinant is %f \n",i,det);  
       cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);  
941      };      };
942        //
943        printf(" done, closing file and exiting\n");
944        //
945    };    };
946        //
   printf(" done, closing file and exiting\n");  
   
947    cf->CloseMatrixFile();    cf->CloseMatrixFile();
948      //
949    cf->Delete();    cf->Delete();
950      //
951  }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23