/[PAMELA software]/calo/flight/FUTILITIES/macros/FCaloADC2MIP.cxx
ViewVC logotype

Diff of /calo/flight/FUTILITIES/macros/FCaloADC2MIP.cxx

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

revision 1.3 by mocchiut, Fri Jun 30 12:08:13 2006 UTC revision 1.5 by mocchiut, Thu Jan 23 11:23:56 2014 UTC
# Line 18  Line 18 
18  #include <sstream>  #include <sstream>
19  #include <string>  #include <string>
20  #include <iostream>  #include <iostream>
21    #include <stdlib.h>
22  #include <TTree.h>  #include <TTree.h>
23  #include <TClassEdit.h>  #include <TClassEdit.h>
24  #include <TObject.h>  #include <TObject.h>
# Line 66  Double_t langaupro(Double_t *); Line 67  Double_t langaupro(Double_t *);
67  #include <FCaloADC2MIPfun.h>  #include <FCaloADC2MIPfun.h>
68  using namespace std;  using namespace std;
69    
70  void FCaloADC2MIP(TString Filename, TString calcalibfile = "", TString OutDir="", TString flist=""){  void FCaloADC2MIP(TString Filename, TString calcalibfile = "", TString OutDir="", TString flist="", Int_t fit=1){
71    const char *pam_calib = pathtocalibration();    const char *pam_calib = pathtocalibration();
72    if ( !strcmp(OutDir.Data(),"") ){    if ( !strcmp(OutDir.Data(),"") ){
73      OutDir = pam_calib;      OutDir = pam_calib;
# Line 486  void FCaloADC2MIP(TString Filename, TStr Line 487  void FCaloADC2MIP(TString Filename, TStr
487      //      //
488      for (Int_t i = 0; i < nevents; i++){      for (Int_t i = 0; i < nevents; i++){
489        evno = i;        evno = i;
490        if ( i%1000 == 0 && i > 0 ) printf(" %iK \n",i/1000);                if ( i%100000 == 0 && i > 0 ) printf(" %i00K \n",i/100000);          
491        //      printf(" %i \n",i);        //      printf(" %i \n",i);
492        //        //
493        // read from the header of the event the absolute time at which it was recorded        // read from the header of the event the absolute time at which it was recorded
# Line 720  void FCaloADC2MIP(TString Filename, TStr Line 721  void FCaloADC2MIP(TString Filename, TStr
721      e++;      e++;
722    };    };
723    //    //
724      if ( !fit ) UPYES=0;
725      //
726    if ( !UPYES ) goto end;    if ( !UPYES ) goto end;
727    printf("\n OK, now I will fit the MIP distribution for each strip \n\n");    printf("\n OK, now I will fit the MIP distribution for each strip \n\n");
728    //    //
# Line 748  void FCaloADC2MIP(TString Filename, TStr Line 751  void FCaloADC2MIP(TString Filename, TStr
751            //            //
752            Double_t fr[2];            Double_t fr[2];
753            Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];            Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
754            if ( SNRPeak > 16. && SNRPeak < 35. ){            //      if ( SNRPeak > 16. && SNRPeak < 35. ){
755              if ( SNRPeak > 16. && SNRPeak < 35. && fp[3] < 8.){
756              fr[0] = (Float_t)SNRPeak - 7.;              fr[0] = (Float_t)SNRPeak - 7.;
757              sv[0]= fp[0];              sv[0]= fp[0];
758              sv[1]= fp[1];              sv[1]= fp[1];
759              sv[2]= fp[2];              sv[2]= fp[2];
760              sv[3]= fp[3];                                sv[3]= fp[3];                  
761            } else {            } else {
762              fr[0] = 19.;              //      fr[0] = 19.;
763                fr[0] = 17.5;
764              sv[0] = 2.8;              sv[0] = 2.8;
765              sv[1] = 21.0;              //      sv[1] = 21.0;
766                sv[1] = 23.0;
767              sv[2] = 1000.0;              sv[2] = 1000.0;
768              sv[3] = 0.1;                                  sv[3] = 0.1;                    
769            };            };
770            Int_t doitagain = 0;            Int_t doitagain = 0;
771          fitting: printf("Fitting strip %i %i %i \n",l,m,n);          fitting: printf("Fitting strip %i %i %i \n",l,m,n);
772            fr[1] = 60.;            //      fr[1] = 60.;
773              fr[1] = 50.;
774            pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.2;            pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.2;
775            plhi[0]=8.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;            //      plhi[0]=8.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;
776              plhi[0]=15.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;
777            Double_t chisqr;            Double_t chisqr;
778            Int_t    ndf;            Int_t    ndf;
779            //            //
# Line 780  void FCaloADC2MIP(TString Filename, TStr Line 788  void FCaloADC2MIP(TString Filename, TStr
788            //            //
789            //            //
790            //            //
791            if ( ( SNRPeak < 0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 3 ){            //      if ( ( SNRPeak < 0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 3 ){
792              if ( ( SNRPeak < 0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 4 ){
793              printf("\n The fit is not good, I will try again, step %i \n\n",doitagain);              printf("\n The fit is not good, I will try again, step %i \n\n",doitagain);
794              doitagain++;              doitagain++;
795              if ( chisqr > 100. ) {              if ( chisqr > 100. ) {
796                fr[0] = fr[0] + 5.;                //              fr[0] = fr[0] + 5.;
797                  fr[0] = fr[0] + 2.;
798                sv[0] = fp[0];                sv[0] = fp[0];
799                sv[1] = fp[1];                sv[1] = fp[1];
800                sv[2] = fp[2];                sv[2] = fp[2];
# Line 2170  void FCalo4224FIT(TString filename = "", Line 2180  void FCalo4224FIT(TString filename = "",
2180    const char *file2;    const char *file2;
2181    file2 = filenome.str().c_str();    file2 = filenome.str().c_str();
2182    TFile *hfile = 0;    TFile *hfile = 0;
2183    hfile = new TFile(file2,"UPDATE","Calorimeter CALIBRATION data");    hfile = new TFile(filevalue.Data(),"UPDATE","Calorimeter CALIBRATION data");
2184    CalorimeterCalibration *calo = 0;    CalorimeterCalibration *calo = 0;
2185    TTree *tree;    TTree *tree;
2186    tree = (TTree*)hfile->Get("CaloADC");      tree = (TTree*)hfile->Get("CaloADC");  
# Line 2184  void FCalo4224FIT(TString filename = "", Line 2194  void FCalo4224FIT(TString filename = "",
2194    file3.str("");    file3.str("");
2195    //  file3 << OutDir.Data() << "/";    //  file3 << OutDir.Data() << "/";
2196    file3 << file;    file3 << file;
2197    TFile hfile3(file3.str().c_str());    TFile hfile3(filename.Data());
2198    Int_t l = 0;    Int_t l = 0;
2199    TCanvas *c1;    TCanvas *c1;
2200    Int_t ty = 1;    Int_t ty = 1;
# Line 2237  void FCalo4224FIT(TString filename = "", Line 2247  void FCalo4224FIT(TString filename = "",
2247          //          //
2248          printf(" **** STRIP: %i %i %i **** \n",l,m,n);          printf(" **** STRIP: %i %i %i **** \n",l,m,n);
2249          printf(" MIP: %f \n",calo->mip[l][m][n]);          printf(" MIP: %f \n",calo->mip[l][m][n]);
2250          printf(" ERMIP: %f \n",calo->ermip[l][m][n]);          //printf(" ERMIP: %f \n",calo->ermip[l][m][n]);
2251          printf(" FP[0-3]: %f / %f / %f / %f \n",calo->fp[0][l][m][n],calo->fp[1][l][m][n],calo->fp[2][l][m][n],calo->fp[3][l][m][n]);          //printf(" FP[0-3]: %f / %f / %f / %f \n",calo->fp[0][l][m][n],calo->fp[1][l][m][n],calo->fp[2][l][m][n],calo->fp[3][l][m][n]);
2252          printf(" FPE[0-3]: %f / %f / %f / %f \n",calo->fpe[0][l][m][n],calo->fpe[1][l][m][n],calo->fpe[2][l][m][n],calo->fpe[3][l][m][n]);          //printf(" FPE[0-3]: %f / %f / %f / %f \n",calo->fpe[0][l][m][n],calo->fpe[1][l][m][n],calo->fpe[2][l][m][n],calo->fpe[3][l][m][n]);
2253          printf(" CHI2: %f \n",calo->chi2[l][m][n]);          //printf(" CHI2: %f \n",calo->chi2[l][m][n]);
2254          printf(" NDF: %f \n",calo->ndf[l][m][n]);          //printf(" NDF: %f \n",calo->ndf[l][m][n]);
2255          printf(" MASK: %f \n",calo->mask[l][m][n]);                      //printf(" MASK: %f \n",calo->mask[l][m][n]);          
2256          //          //
2257          // what to do?          // what to do?
2258          //          //
# Line 2312  void FCalo4224FIT(TString filename = "", Line 2322  void FCalo4224FIT(TString filename = "",
2322                input.str("");                input.str("");
2323                input << tellme;                input << tellme;
2324                out.str("y");                out.str("y");
2325                if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {                                  if ( !strcmp(input.str().c_str(),out.str().c_str()) ) {                  
2326                  calo->mip[l][m][n] = (float)SNRPeak;                  calo->mip[l][m][n] = (float)SNRPeak;
2327                  calo->ermip[l][m][n] = (float)ffpe[1];                  calo->ermip[l][m][n] = (float)ffpe[1];
2328                  for (Int_t a = 0; a < 4 ; a++){                  for (Int_t a = 0; a < 4 ; a++){
# Line 2454  void FCalo4224MIPVALUES(TString filevalu Line 2464  void FCalo4224MIPVALUES(TString filevalu
2464    const char *file2;    const char *file2;
2465    file2 = filenome.str().c_str();    file2 = filenome.str().c_str();
2466    TFile *hfile = 0;    TFile *hfile = 0;
2467    hfile = new TFile(file2,"READ","Calorimeter CALIBRATION data");    hfile = new TFile(filevalue.Data(),"READ","Calorimeter CALIBRATION data");
2468    CalorimeterCalibration *calo = 0;    CalorimeterCalibration *calo = 0;
2469    TTree *tree;    TTree *tree;
2470    tree = (TTree*)hfile->Get("CaloADC");      tree = (TTree*)hfile->Get("CaloADC");  
# Line 2498  void FCalo4224MIPVALUES(TString filevalu Line 2508  void FCalo4224MIPVALUES(TString filevalu
2508          mg->Add(graph);          mg->Add(graph);
2509          printf(" **** STRIP: %i %i %i **** \n",l,m,n);          printf(" **** STRIP: %i %i %i **** \n",l,m,n);
2510          printf(" MIP: %f \n",calo->mip[l][m][n]);          printf(" MIP: %f \n",calo->mip[l][m][n]);
2511          printf(" ERMIP: %f \n",calo->ermip[l][m][n]);          //printf(" ERMIP: %f \n",calo->ermip[l][m][n]);
2512          printf(" FP[0-3]: %f / %f / %f / %f \n",calo->fp[0][l][m][n],calo->fp[1][l][m][n],calo->fp[2][l][m][n],calo->fp[3][l][m][n]);          //printf(" FP[0-3]: %f / %f / %f / %f \n",calo->fp[0][l][m][n],calo->fp[1][l][m][n],calo->fp[2][l][m][n],calo->fp[3][l][m][n]);
2513          printf(" FPE[0-3]: %f / %f / %f / %f \n",calo->fpe[0][l][m][n],calo->fpe[1][l][m][n],calo->fpe[2][l][m][n],calo->fpe[3][l][m][n]);          //printf(" FPE[0-3]: %f / %f / %f / %f \n",calo->fpe[0][l][m][n],calo->fpe[1][l][m][n],calo->fpe[2][l][m][n],calo->fpe[3][l][m][n]);
2514          printf(" CHI2: %f \n",calo->chi2[l][m][n]);          //printf(" CHI2: %f \n",calo->chi2[l][m][n]);
2515          printf(" NDF: %f \n",calo->ndf[l][m][n]);          //printf(" NDF: %f \n",calo->ndf[l][m][n]);
2516          printf(" MASK: %f \n",calo->mask[l][m][n]);                      //printf(" MASK: %f \n",calo->mask[l][m][n]);          
2517          n++;          n++;
2518        };        };
2519        l++;        l++;
# Line 2841  void FCaloBAKFIT(TString filename, TStri Line 2851  void FCaloBAKFIT(TString filename, TStri
2851          bmippa << " " << n;                      bmippa << " " << n;            
2852          pfmip[l][m][n] = (TH1F*)temp->Clone(bmippa.str().c_str());          pfmip[l][m][n] = (TH1F*)temp->Clone(bmippa.str().c_str());
2853          Int_t doitagain = 0;          Int_t doitagain = 0;
2854            Int_t ndoitagain = 0;
2855          Double_t fr[2];          Double_t fr[2];
2856          Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];          Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
2857          if ( SNRPeak > 16. && SNRPeak < 35. ){          //      if ( SNRPeak > 16. && SNRPeak < 35. ){
2858            if ( SNRPeak > 16. && SNRPeak < 35. && fp[3] < 8. ){
2859            fr[0] = (Float_t)SNRPeak - 7.;            fr[0] = (Float_t)SNRPeak - 7.;
2860            sv[0]= fp[0];            sv[0]= fp[0];
2861            sv[1]= fp[1];            sv[1]= fp[1];
2862            sv[2]= fp[2];            sv[2]= fp[2];
2863            sv[3]= fp[3];                        sv[3]= fp[3];            
2864          } else {          } else {
2865            fr[0] = 19.;            fr[0] = 12.5;
2866            sv[0] = 2.8;            sv[0] = 2.8;
2867            sv[1] = 21.0;            sv[1] = 23.0;
2868            sv[2] = 1000.0;            sv[2] = 1000.0;
2869            sv[3] = 0.1;                          sv[3] = 0.1;              
2870          };          };
2871            fr[1] = 50.;
2872            //      fr[1] = 60.;
2873        fitting: printf("Fitting strip %i %i %i \n",l,m,n);        fitting: printf("Fitting strip %i %i %i \n",l,m,n);
         fr[1] = 60.;  
2874          pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.2;          pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.2;
2875          plhi[0]=8.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;          plhi[0]=15.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;
2876          Double_t chisqr;          Double_t chisqr;
2877          Int_t    ndf;          Int_t    ndf;
2878          pfitsnr[l][m][n] = langaufit(pfmip[l][m][n],fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf,"QR0");          pfitsnr[l][m][n] = langaufit(pfmip[l][m][n],fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf,"QR0");
# Line 2867  void FCaloBAKFIT(TString filename, TStri Line 2880  void FCaloBAKFIT(TString filename, TStri
2880          Double_t SNRPeak = langaupro(fp);          Double_t SNRPeak = langaupro(fp);
2881          printf("\n Conversion factor: %f \n\n",SNRPeak);          printf("\n Conversion factor: %f \n\n",SNRPeak);
2882          //          //
2883            c1->cd();                  
2884            pfmip[l][m][n]->DrawCopy();        
2885            pfitsnr[l][m][n]->DrawCopy("lsame");
2886            c1->Modified();
2887            c1->Update();
2888          //          //
2889          //          //
2890          if ( (SNRPeak<0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 3 ){          if ( (SNRPeak<0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 11 ){
2891            printf("\n The fit is not good, I will try again, step %i \n\n",doitagain);            printf("\n The fit is not good, I will try again, step %i \n\n",doitagain);
2892            doitagain++;            doitagain++;
2893            if ( chisqr > 100. ) {            if ( chisqr > 100. ) {
2894              fr[0] = fr[0] + 5.;              if ( fr[0] < 20. && !(doitagain%2) ) fr[0] += 1.;
2895                if ( fr[1] > 35. && doitagain%2 ) fr[1] -= 3.;
2896                //      if ( doitagain > 1 ) fr[1] = 35.;
2897              sv[0] = fp[0];              sv[0] = fp[0];
2898              sv[1] = fp[1];              sv[1] = fp[1];
2899              sv[2] = fp[2];              sv[2] = fp[2];
2900              sv[3] = fp[3];                                sv[3] = fp[3];      
2901                printf(" chi \n");
2902              goto fitting;              goto fitting;
2903            };            };
2904            if ( doitagain == 1 ) {            ndoitagain++;
2905              if ( ndoitagain == 1 ) {
2906              fr[0] = 19.;              fr[0] = 19.;
2907              sv[0] = 2.8;              sv[0] = 2.8;
2908              sv[1] = 21.0;              sv[1] = 21.0;
2909              sv[2] = 1000.0;              sv[2] = 1000.0;
2910              sv[3] = 0.1;              sv[3] = 0.1;
2911                printf(" uno \n");
2912              goto fitting;              goto fitting;
2913            };            };
2914            if ( doitagain == 2 ) {            if ( ndoitagain == 2 ) {
2915              fr[0] = 22.;              fr[0] = 22.;
2916              sv[0] = 2.8;              sv[0] = 2.8;
2917              sv[1] = 25.0;              sv[1] = 25.0;
2918              sv[2] = 1000.0;              sv[2] = 1000.0;
2919              sv[3] = 0.1;              sv[3] = 0.1;
2920                printf(" due \n");
2921              goto fitting;              goto fitting;
2922            };            };
2923            if ( doitagain == 3 ) {            if ( ndoitagain == 3 ) {
2924              fr[0] = 12.;              fr[0] = 12.;
2925              sv[0] = 2.8;              sv[0] = 2.8;
2926              sv[1] = 15.0;              sv[1] = 15.0;
2927              sv[2] = 1000.0;              sv[2] = 1000.0;
2928              sv[3] = 0.1;              sv[3] = 0.1;
2929                printf(" tre \n");
2930              goto fitting;              goto fitting;
2931            };            };
2932          };          };
# Line 3000  void FCaloROOT2TXT(TString rootple, TStr Line 3025  void FCaloROOT2TXT(TString rootple, TStr
3025      for (Int_t k = 0; k < 22; k++ ){      for (Int_t k = 0; k < 22; k++ ){
3026        for (Int_t l = 0; l < 96; l++ ){        for (Int_t l = 0; l < 96; l++ ){
3027          mip = 0.;          mip = 0.;
3028          if ( (ccalo->fp[1][m][k][l] > 20. && ccalo->fp[1][m][k][l] < 32.) || ccalo->mask[m][k][l] == 1. ) {          //      if ( (ccalo->fp[1][m][k][l] > 20. && ccalo->fp[1][m][k][l] < 32.) || ccalo->mask[m][k][l] == 1. ) {
3029    //      if ( (ccalo->fp[1][m][k][l] > 15. && ccalo->fp[1][m][k][l] < 45.) || ccalo->mask[m][k][l] == 1. ) {
3030            mip = ccalo->mip[m][k][l];            mip = ccalo->mip[m][k][l];
3031          } else {  //      } else {
3032            mip = 26. ;            //      mip = 26. ;
3033          };  //        mip = 28. ;
3034          if ( mip == 0 ) mip = 26. ;  //      };
3035            //      if ( mip == 0 ) mip = 26. ;
3036    //      if ( mip == 0 ) mip = 28. ;
3037          //      if ( m == 1 && k == 17 ) printf(" %i mip = %f \n",l,mip);          //      if ( m == 1 && k == 17 ) printf(" %i mip = %f \n",l,mip);
3038          fwrite(&mip,sizeof(mip),1,f);          fwrite(&mip,sizeof(mip),1,f);
3039        };        };
# Line 3018  void FCaloREADTXT(TString txtple){ Line 3046  void FCaloREADTXT(TString txtple){
3046    FILE *f;    FILE *f;
3047    f = fopen(txtple.Data(),"rb");    f = fopen(txtple.Data(),"rb");
3048    Float_t mip = 0.;    Float_t mip = 0.;
3049      Float_t tmip = 0.;
3050    for (Int_t m = 0; m < 2 ; m++ ){    for (Int_t m = 0; m < 2 ; m++ ){
3051      for (Int_t k = 0; k < 22; k++ ){      for (Int_t k = 0; k < 22; k++ ){
3052        for (Int_t l = 0; l < 96; l++ ){        for (Int_t l = 0; l < 96; l++ ){
3053          mip = 0.;          mip = 0.;
3054          fread(&mip,sizeof(mip),1,f);          fread(&mip,sizeof(mip),1,f);
3055            tmip += mip;
3056            printf(" m %i k %i l %i mip %f \n",m,k,l,mip);
3057          //      if ( m == 1 && k == 17 ) printf(" %i mip = %f \n",l,mip);          //      if ( m == 1 && k == 17 ) printf(" %i mip = %f \n",l,mip);
3058        };        };
3059      };      };
3060    };    };
3061    fclose(f);    fclose(f);
3062      //
3063      printf(" Average mip is %f \n",tmip/4224.);
3064      //
3065  }  }
3066    
3067  void FCaloALIG2DAT(TString txtple){  void FCaloALIG2DAT(TString txtple){
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23