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

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

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

revision 1.3 by mocchiut, Thu Dec 13 17:08:11 2007 UTC revision 1.4 by mocchiut, Tue Dec 18 09:55:07 2007 UTC
# Line 18  CaloFranzini::CaloFranzini(PamLevel2 *l2 Line 18  CaloFranzini::CaloFranzini(PamLevel2 *l2
18    //    //
19    file = NULL;    file = NULL;
20    brig = NULL;    brig = NULL;
21      brigm = NULL;
22    nbin = 0;    nbin = 0;
23    //    //
24    L2 = l2p;    L2 = l2p;
# Line 209  void CaloFranzini::Process(Int_t itr){ Line 210  void CaloFranzini::Process(Int_t itr){
210      degfre = TMath::Min(dgf,NC);      degfre = TMath::Min(dgf,NC);
211      //      //
212      if ( degfre > 0 ){      if ( degfre > 0 ){
       TArrayF *qplmean = this->LoadLongAverage(rig);  
       TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig);  
213        for (Int_t i = 0; i < degfre; i++){        for (Int_t i = 0; i < degfre; i++){
214          if ( i != mask18b ){          if ( i != mask18b ){
215            for (Int_t j = 0; j < degfre; j++){              for (Int_t j = 0; j < degfre; j++){  
216              if ( j != mask18b ){              if ( j != mask18b ){
217                longtzeta += (qplane[i]-qplmean->At(i)) * (*lmatrix)[i][j] * (qplane[j]-qplmean->At(j));                longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
218              };              };
219            };            };
220          };          };
# Line 262  Bool_t CaloFranzini::CreateMatrixFile(TS Line 261  Bool_t CaloFranzini::CreateMatrixFile(TS
261    //    //
262  }  }
263    
264    Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){  
265      //
266      file = new TFile(matrixfile.Data(),"UPDATE");
267      //
268      if ( !file || file->IsZombie() ){
269        printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
270        return(false);
271      };
272      //
273      return(true);
274      //
275    }
276    
277  void CaloFranzini::WriteNumBin(Int_t numbin){  void CaloFranzini::WriteNumBin(Int_t numbin){
278    file->cd();    file->cd();
279    TArrayI nbi(1, &numbin);    TArrayI nbi(1, &numbin);
# Line 316  Bool_t CaloFranzini::Open(TString matrix Line 328  Bool_t CaloFranzini::Open(TString matrix
328      return(false);      return(false);
329    };    };
330    //    //
331      if ( !this->LoadBin() ){
332        printf(" %s \n",matrixfile.Data());
333        return(false);
334      };
335      if ( !this->LoadMatrices() ){
336        printf(" %s \n",matrixfile.Data());
337        return(false);
338      };
339      //
340      //
341      return(true);
342      //
343    }
344    
345    Bool_t CaloFranzini::LoadBin(){
346      //
347    TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");    TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");
348    if ( !numbin ){    if ( !numbin ){
349      printf(" ERROR: cannot read number of bins from file %s \n",matrixfile.Data());      printf(" ERROR: cannot read number of bins from file ");
350      return(false);      return(false);
351    };    };
352    nbin = (Int_t)numbin->At(0);    nbin = (Int_t)numbin->At(0);
353    if ( nbin <= 0 ){    if ( nbin <= 0 ){
354      printf(" ERROR: cannot work with 0 energy bins (from file %s) \n",matrixfile.Data());      printf(" ERROR: cannot work with 0 energy bins from file ");
355      return(false);      return(false);
356    };    };
357    //    //
358    brig = (TArrayF*)file->Get("binrig");    brig = (TArrayF*)file->Get("binrig");
359    if ( !brig ){    if ( !brig ){
360      printf(" ERROR: cannot read rigidity binning from file %s \n",matrixfile.Data());      printf(" ERROR: cannot read rigidity binning from file ");
361      return(false);      return(false);
362    };    };
363    //    //
364      brigm=(TArrayF*)file->Get("binrigmean");
365      if ( !brigm ){
366        printf(" ERROR: cannot read mean rigidity binning from file ");
367        return(false);
368      };
369      //
370      for (Int_t i=0;i<17;i++){
371            TString name = Form("qplmeann%i",i);
372            qplmean[i] = (TArrayF*)file->Get(name.Data());
373            if ( !qplmean[i] ){
374                    printf(" ERROR: cannot read average from file ");
375                    return(false);
376            };
377      };
378      //
379    return(true);    return(true);
380    }
381    
382    Bool_t CaloFranzini::LoadMatrices(){
383    //    //
384      for (Int_t i=0;i<17;i++){
385            TString name1 = Form("matrixn%i",i);
386            hmat[i] = (TMatrixD*)file->Get(name1.Data());
387      };
388      //
389      return(true);
390  }  }
391    
392  TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){  TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
393    //    //
394    TString name;    Int_t mv = 0;
395    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
396      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
397        name = Form("matrixn%i",i);        mv = i;
398        break;        break;
399      };      };
400    };    };
401    if ( rig < brig->At(0) ){    if ( rig < brig->At(0) ){
402      printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));      printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
403      name = "matrixn0";      mv = 0;
     printf(" Using matrix %s \n",name.Data());  
404    };    };
405    if ( rig >= brig->At(nbin-1) ){    if ( rig >= brig->At(nbin-1) ){
406      printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));      printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
407      name = Form("matrixn%i",nbin-2);      mv = nbin-2;
     printf(" Using matrix %s \n",name.Data());  
408    };    };
409    //    //
410    TMatrixD *matrix = (TMatrixD*)file->Get(name.Data());    return(hmat[mv]);
   //  
   return(matrix);  
411    //    //
412  }  }
413    
414    
415  TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){  TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
416    //    //
417    TString name;    Int_t mv=0;
418    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
419      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
420        name = Form("qplmeann%i",i);          mv = i;
421        break;        break;
422      };      };
423    };    };
424    if ( rig < brig->At(0) ){    if ( rig < brig->At(0) ){
425      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
426      name = "qplmeann0";      mv = 0;
     printf(" Using qplmean %s \n",name.Data());  
427    };    };
428    if ( rig >= brig->At(nbin-1) ){    if ( rig >= brig->At(nbin-1) ){
429      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
430      name = Form("qplmeann%i",nbin-2);      mv=nbin-2;
     printf(" Using qplmean %s \n",name.Data());  
431    };    };
432    //    //
433    TArrayF *qplmean = (TArrayF*)file->Get(name.Data());    return(qplmean[mv]);
   //  
   return(qplmean);  
434    //    //
435  }  }
436    
437    Float_t CaloFranzini::GetAverageAt(Int_t plane, Float_t rig){
438      //
439      Int_t therigb = 0;
440      for (Int_t i = 0; i<nbin-2; i++){
441        if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
442          therigb = i;
443          break;
444        };
445      };
446      //
447      Float_t minrig;
448      Float_t maxrig;
449      //
450      //
451      if ( rig < brigm->At(0) ){
452        if ( rig < brig->At(0) ){
453          printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
454        };
455        therigb = 0;
456      };
457      if ( rig >= brigm->At(nbin-2) ){
458        if ( rig >= brig->At(nbin-2) ) {
459          printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
460        };
461        therigb = nbin - 3;
462      };
463      //
464      minrig = brigm->At(therigb);
465      maxrig = brigm->At(therigb+1);
466      //
467      Float_t minene = (*qplmean[therigb])[plane];
468      Float_t maxene = (*qplmean[therigb+1])[plane];
469      //
470      if ( maxrig == minrig ){
471       printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
472       return(0.);
473      };
474      Float_t b = log(maxene/minene)/(maxrig-minrig);
475      Float_t a = minene/exp(minrig*b);
476      Float_t ave = a*exp(b*rig);
477      //
478      return(ave);
479      //
480    }
481    
482    Float_t CaloFranzini::GetHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
483      Int_t therigb = 0;
484      for (Int_t i = 0; i<nbin-2; i++){
485        if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
486          therigb = i;
487          break;
488        };
489      };
490      //
491      Float_t minrig;
492      Float_t maxrig;
493      //
494      if ( rig < brigm->At(0) ){
495        if ( rig < brig->At(0) ){
496          printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
497        };
498        therigb = 0;
499      };
500      if ( rig >= brigm->At(nbin-2) ){
501        if ( rig >= brig->At(nbin-2) ) {
502          printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
503        };
504        therigb = nbin - 3;
505      };
506      //
507      minrig = brigm->At(therigb);
508      maxrig = brigm->At(therigb+1);
509      //
510      Float_t minene = (*hmat[therigb])[iindex][jindex];
511      Float_t maxene = (*hmat[therigb+1])[iindex][jindex];
512      //
513      Float_t b = log(maxene/minene)/(maxrig-minrig);
514      Float_t a = minene/exp(minrig*b);
515      Float_t ave = a*exp(b*rig);
516      //
517      return(ave);
518      //
519    }
520    
521  void CaloFranzini::DrawLongAverage(Float_t rig){  void CaloFranzini::DrawLongAverage(Float_t rig){
522    //    //
# Line 412  void CaloFranzini::DrawLongAverage(Float Line 538  void CaloFranzini::DrawLongAverage(Float
538    if ( thf ) thf->Delete();    if ( thf ) thf->Delete();
539    thf = new TH1F(thid,thid,44,-0.5,43.5);    thf = new TH1F(thid,thid,44,-0.5,43.5);
540    tcf->cd();    tcf->cd();
   //  Int_t pp=0;  
541    Float_t qpl[43];    Float_t qpl[43];
542    for (Int_t st=0;st<43;st++){    for (Int_t st=0;st<43;st++){
543      qpl[st] = ll->At(st);      qpl[st] = ll->At(st);
544        printf("st %i qpl %f\n",st,qpl[st]);
545    };    };
546    for (Int_t st=0;st<44;st++){    for (Int_t st=0;st<44;st++){
547      if ( st == 37 ){      if ( st == 37 ){
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.23