/[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.1 by mocchiut, Wed Mar 22 15:07:47 2006 UTC revision 1.4 by mocchiut, Mon Sep 22 20:18:43 2008 UTC
# Line 1  Line 1 
1  //  //
2  //- Emiliano Mocchiutti  //- Emiliano Mocchiutti
3  //  //
4  //   FCaloADC2MIP.c      version 1.01  (2006-03-20)  //   FCaloADC2MIP.c      version 1.02  (2006-03-31)
5  //  //
6  //   The only input needed is  //   The only input needed is
7  //  //
8  //   Changelog:  //   Changelog:
9  //  //
10    //   1.01 - 1.02 (2006-03-31): Added FCaloALIG2DAT and FCaloREADALIG to create and read the calorimeter alignment data file.
11    //
12  //   1.00 - 1.01 (2006-03-20): First flight version, read single yoda file.  //   1.00 - 1.01 (2006-03-20): First flight version, read single yoda file.
13  //  //
14  //   0.00 - 1.00 (2006-03-20): Clone of CaloADC2MIP v4r01.  //   0.00 - 1.00 (2006-03-20): Clone of CaloADC2MIP v4r01.
# Line 64  Double_t langaupro(Double_t *); Line 66  Double_t langaupro(Double_t *);
66  #include <FCaloADC2MIPfun.h>  #include <FCaloADC2MIPfun.h>
67  using namespace std;  using namespace std;
68    
69  void FCaloADC2MIP(TString Filename, TString calcalibfile = "", TString OutDir="", TString flist=""){  void FCaloADC2MIP(TString Filename, TString calcalibfile = "", TString OutDir="", TString flist="", Int_t fit=1){
70    const char *pam_calib = pathtocalibration();    const char *pam_calib = pathtocalibration();
71    if ( !strcmp(OutDir.Data(),"") ){    if ( !strcmp(OutDir.Data(),"") ){
72      OutDir = pam_calib;      OutDir = pam_calib;
# Line 484  void FCaloADC2MIP(TString Filename, TStr Line 486  void FCaloADC2MIP(TString Filename, TStr
486      //      //
487      for (Int_t i = 0; i < nevents; i++){      for (Int_t i = 0; i < nevents; i++){
488        evno = i;        evno = i;
489        if ( i%1000 == 0 && i > 0 ) printf(" %iK \n",i/1000);                if ( i%100000 == 0 && i > 0 ) printf(" %i00K \n",i/100000);          
490        //      printf(" %i \n",i);        //      printf(" %i \n",i);
491        //        //
492        // 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 718  void FCaloADC2MIP(TString Filename, TStr Line 720  void FCaloADC2MIP(TString Filename, TStr
720      e++;      e++;
721    };    };
722    //    //
723      if ( !fit ) UPYES=0;
724      //
725    if ( !UPYES ) goto end;    if ( !UPYES ) goto end;
726    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");
727    //    //
# Line 746  void FCaloADC2MIP(TString Filename, TStr Line 750  void FCaloADC2MIP(TString Filename, TStr
750            //            //
751            Double_t fr[2];            Double_t fr[2];
752            Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];            Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
753            if ( SNRPeak > 16. && SNRPeak < 35. ){            //      if ( SNRPeak > 16. && SNRPeak < 35. ){
754              if ( SNRPeak > 16. && SNRPeak < 35. && fp[3] < 8.){
755              fr[0] = (Float_t)SNRPeak - 7.;              fr[0] = (Float_t)SNRPeak - 7.;
756              sv[0]= fp[0];              sv[0]= fp[0];
757              sv[1]= fp[1];              sv[1]= fp[1];
758              sv[2]= fp[2];              sv[2]= fp[2];
759              sv[3]= fp[3];                                sv[3]= fp[3];                  
760            } else {            } else {
761              fr[0] = 19.;              //      fr[0] = 19.;
762                fr[0] = 17.5;
763              sv[0] = 2.8;              sv[0] = 2.8;
764              sv[1] = 21.0;              //      sv[1] = 21.0;
765                sv[1] = 23.0;
766              sv[2] = 1000.0;              sv[2] = 1000.0;
767              sv[3] = 0.1;                                  sv[3] = 0.1;                    
768            };            };
769            Int_t doitagain = 0;            Int_t doitagain = 0;
770          fitting: printf("Fitting strip %i %i %i \n",l,m,n);          fitting: printf("Fitting strip %i %i %i \n",l,m,n);
771            fr[1] = 60.;            //      fr[1] = 60.;
772              fr[1] = 50.;
773            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;
774            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;
775              plhi[0]=15.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;
776            Double_t chisqr;            Double_t chisqr;
777            Int_t    ndf;            Int_t    ndf;
778            //            //
# Line 778  void FCaloADC2MIP(TString Filename, TStr Line 787  void FCaloADC2MIP(TString Filename, TStr
787            //            //
788            //            //
789            //            //
790            if ( ( SNRPeak < 0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 3 ){            //      if ( ( SNRPeak < 0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 3 ){
791              if ( ( SNRPeak < 0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 4 ){
792              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);
793              doitagain++;              doitagain++;
794              if ( chisqr > 100. ) {              if ( chisqr > 100. ) {
795                fr[0] = fr[0] + 5.;                //              fr[0] = fr[0] + 5.;
796                  fr[0] = fr[0] + 2.;
797                sv[0] = fp[0];                sv[0] = fp[0];
798                sv[1] = fp[1];                sv[1] = fp[1];
799                sv[2] = fp[2];                sv[2] = fp[2];
# Line 2168  void FCalo4224FIT(TString filename = "", Line 2179  void FCalo4224FIT(TString filename = "",
2179    const char *file2;    const char *file2;
2180    file2 = filenome.str().c_str();    file2 = filenome.str().c_str();
2181    TFile *hfile = 0;    TFile *hfile = 0;
2182    hfile = new TFile(file2,"UPDATE","Calorimeter CALIBRATION data");    hfile = new TFile(filevalue.Data(),"UPDATE","Calorimeter CALIBRATION data");
2183    CalorimeterCalibration *calo = 0;    CalorimeterCalibration *calo = 0;
2184    TTree *tree;    TTree *tree;
2185    tree = (TTree*)hfile->Get("CaloADC");      tree = (TTree*)hfile->Get("CaloADC");  
# Line 2182  void FCalo4224FIT(TString filename = "", Line 2193  void FCalo4224FIT(TString filename = "",
2193    file3.str("");    file3.str("");
2194    //  file3 << OutDir.Data() << "/";    //  file3 << OutDir.Data() << "/";
2195    file3 << file;    file3 << file;
2196    TFile hfile3(file3.str().c_str());    TFile hfile3(filename.Data());
2197    Int_t l = 0;    Int_t l = 0;
2198    TCanvas *c1;    TCanvas *c1;
2199    Int_t ty = 1;    Int_t ty = 1;
# Line 2310  void FCalo4224FIT(TString filename = "", Line 2321  void FCalo4224FIT(TString filename = "",
2321                input.str("");                input.str("");
2322                input << tellme;                input << tellme;
2323                out.str("y");                out.str("y");
2324                if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {                                  if ( !strcmp(input.str().c_str(),out.str().c_str()) ) {                  
2325                  calo->mip[l][m][n] = (float)SNRPeak;                  calo->mip[l][m][n] = (float)SNRPeak;
2326                  calo->ermip[l][m][n] = (float)ffpe[1];                  calo->ermip[l][m][n] = (float)ffpe[1];
2327                  for (Int_t a = 0; a < 4 ; a++){                  for (Int_t a = 0; a < 4 ; a++){
# Line 2452  void FCalo4224MIPVALUES(TString filevalu Line 2463  void FCalo4224MIPVALUES(TString filevalu
2463    const char *file2;    const char *file2;
2464    file2 = filenome.str().c_str();    file2 = filenome.str().c_str();
2465    TFile *hfile = 0;    TFile *hfile = 0;
2466    hfile = new TFile(file2,"READ","Calorimeter CALIBRATION data");    hfile = new TFile(filevalue.Data(),"READ","Calorimeter CALIBRATION data");
2467    CalorimeterCalibration *calo = 0;    CalorimeterCalibration *calo = 0;
2468    TTree *tree;    TTree *tree;
2469    tree = (TTree*)hfile->Get("CaloADC");      tree = (TTree*)hfile->Get("CaloADC");  
# Line 2839  void FCaloBAKFIT(TString filename, TStri Line 2850  void FCaloBAKFIT(TString filename, TStri
2850          bmippa << " " << n;                      bmippa << " " << n;            
2851          pfmip[l][m][n] = (TH1F*)temp->Clone(bmippa.str().c_str());          pfmip[l][m][n] = (TH1F*)temp->Clone(bmippa.str().c_str());
2852          Int_t doitagain = 0;          Int_t doitagain = 0;
2853            Int_t ndoitagain = 0;
2854          Double_t fr[2];          Double_t fr[2];
2855          Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];          Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
2856          if ( SNRPeak > 16. && SNRPeak < 35. ){          //      if ( SNRPeak > 16. && SNRPeak < 35. ){
2857            if ( SNRPeak > 16. && SNRPeak < 35. && fp[3] < 8. ){
2858            fr[0] = (Float_t)SNRPeak - 7.;            fr[0] = (Float_t)SNRPeak - 7.;
2859            sv[0]= fp[0];            sv[0]= fp[0];
2860            sv[1]= fp[1];            sv[1]= fp[1];
2861            sv[2]= fp[2];            sv[2]= fp[2];
2862            sv[3]= fp[3];                        sv[3]= fp[3];            
2863          } else {          } else {
2864            fr[0] = 19.;            fr[0] = 12.5;
2865            sv[0] = 2.8;            sv[0] = 2.8;
2866            sv[1] = 21.0;            sv[1] = 23.0;
2867            sv[2] = 1000.0;            sv[2] = 1000.0;
2868            sv[3] = 0.1;                          sv[3] = 0.1;              
2869          };          };
2870            fr[1] = 50.;
2871            //      fr[1] = 60.;
2872        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.;  
2873          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;
2874          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;
2875          Double_t chisqr;          Double_t chisqr;
2876          Int_t    ndf;          Int_t    ndf;
2877          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 2865  void FCaloBAKFIT(TString filename, TStri Line 2879  void FCaloBAKFIT(TString filename, TStri
2879          Double_t SNRPeak = langaupro(fp);          Double_t SNRPeak = langaupro(fp);
2880          printf("\n Conversion factor: %f \n\n",SNRPeak);          printf("\n Conversion factor: %f \n\n",SNRPeak);
2881          //          //
2882            c1->cd();                  
2883            pfmip[l][m][n]->DrawCopy();        
2884            pfitsnr[l][m][n]->DrawCopy("lsame");
2885            c1->Modified();
2886            c1->Update();
2887          //          //
2888          //          //
2889          if ( (SNRPeak<0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 3 ){          if ( (SNRPeak<0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 11 ){
2890            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);
2891            doitagain++;            doitagain++;
2892            if ( chisqr > 100. ) {            if ( chisqr > 100. ) {
2893              fr[0] = fr[0] + 5.;              if ( fr[0] < 20. && !(doitagain%2) ) fr[0] += 1.;
2894                if ( fr[1] > 35. && doitagain%2 ) fr[1] -= 3.;
2895                //      if ( doitagain > 1 ) fr[1] = 35.;
2896              sv[0] = fp[0];              sv[0] = fp[0];
2897              sv[1] = fp[1];              sv[1] = fp[1];
2898              sv[2] = fp[2];              sv[2] = fp[2];
2899              sv[3] = fp[3];                                sv[3] = fp[3];      
2900                printf(" chi \n");
2901              goto fitting;              goto fitting;
2902            };            };
2903            if ( doitagain == 1 ) {            ndoitagain++;
2904              if ( ndoitagain == 1 ) {
2905              fr[0] = 19.;              fr[0] = 19.;
2906              sv[0] = 2.8;              sv[0] = 2.8;
2907              sv[1] = 21.0;              sv[1] = 21.0;
2908              sv[2] = 1000.0;              sv[2] = 1000.0;
2909              sv[3] = 0.1;              sv[3] = 0.1;
2910                printf(" uno \n");
2911              goto fitting;              goto fitting;
2912            };            };
2913            if ( doitagain == 2 ) {            if ( ndoitagain == 2 ) {
2914              fr[0] = 22.;              fr[0] = 22.;
2915              sv[0] = 2.8;              sv[0] = 2.8;
2916              sv[1] = 25.0;              sv[1] = 25.0;
2917              sv[2] = 1000.0;              sv[2] = 1000.0;
2918              sv[3] = 0.1;              sv[3] = 0.1;
2919                printf(" due \n");
2920              goto fitting;              goto fitting;
2921            };            };
2922            if ( doitagain == 3 ) {            if ( ndoitagain == 3 ) {
2923              fr[0] = 12.;              fr[0] = 12.;
2924              sv[0] = 2.8;              sv[0] = 2.8;
2925              sv[1] = 15.0;              sv[1] = 15.0;
2926              sv[2] = 1000.0;              sv[2] = 1000.0;
2927              sv[3] = 0.1;              sv[3] = 0.1;
2928                printf(" tre \n");
2929              goto fitting;              goto fitting;
2930            };            };
2931          };          };
# Line 2998  void FCaloROOT2TXT(TString rootple, TStr Line 3024  void FCaloROOT2TXT(TString rootple, TStr
3024      for (Int_t k = 0; k < 22; k++ ){      for (Int_t k = 0; k < 22; k++ ){
3025        for (Int_t l = 0; l < 96; l++ ){        for (Int_t l = 0; l < 96; l++ ){
3026          mip = 0.;          mip = 0.;
3027          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. ) {
3028    //      if ( (ccalo->fp[1][m][k][l] > 15. && ccalo->fp[1][m][k][l] < 45.) || ccalo->mask[m][k][l] == 1. ) {
3029            mip = ccalo->mip[m][k][l];            mip = ccalo->mip[m][k][l];
3030          } else {  //      } else {
3031            mip = 26. ;            //      mip = 26. ;
3032          };  //        mip = 28. ;
3033          if ( mip == 0 ) mip = 26. ;  //      };
3034            //      if ( mip == 0 ) mip = 26. ;
3035    //      if ( mip == 0 ) mip = 28. ;
3036          //      if ( m == 1 && k == 17 ) printf(" %i mip = %f \n",l,mip);          //      if ( m == 1 && k == 17 ) printf(" %i mip = %f \n",l,mip);
3037          fwrite(&mip,sizeof(mip),1,f);          fwrite(&mip,sizeof(mip),1,f);
3038        };        };
# Line 3016  void FCaloREADTXT(TString txtple){ Line 3045  void FCaloREADTXT(TString txtple){
3045    FILE *f;    FILE *f;
3046    f = fopen(txtple.Data(),"rb");    f = fopen(txtple.Data(),"rb");
3047    Float_t mip = 0.;    Float_t mip = 0.;
3048      Float_t tmip = 0.;
3049    for (Int_t m = 0; m < 2 ; m++ ){    for (Int_t m = 0; m < 2 ; m++ ){
3050      for (Int_t k = 0; k < 22; k++ ){      for (Int_t k = 0; k < 22; k++ ){
3051        for (Int_t l = 0; l < 96; l++ ){        for (Int_t l = 0; l < 96; l++ ){
3052          mip = 0.;          mip = 0.;
3053          fread(&mip,sizeof(mip),1,f);          fread(&mip,sizeof(mip),1,f);
3054            tmip += mip;
3055            printf(" m %i k %i l %i mip %f \n",m,k,l,mip);
3056          //      if ( m == 1 && k == 17 ) printf(" %i mip = %f \n",l,mip);          //      if ( m == 1 && k == 17 ) printf(" %i mip = %f \n",l,mip);
3057        };        };
3058      };      };
3059    };    };
3060    fclose(f);    fclose(f);
3061      //
3062      printf(" Average mip is %f \n",tmip/4224.);
3063      //
3064    }
3065    
3066    void FCaloALIG2DAT(TString txtple){
3067    //   clevel1_.xalig = 121.1;
3068    //   clevel1_.yalig = 119.8;
3069    //   clevel1_.zalig = -263.1;
3070      char value[256];
3071      Int_t ixalig;
3072      Int_t iyalig;
3073      Int_t izalig;
3074      Int_t iemip;
3075      Float_t xalig;
3076      Float_t yalig;
3077      Float_t zalig;
3078      Float_t emip;
3079      //
3080      FILE *f;
3081      f = fopen(txtple.Data(),"wb");
3082      //
3083      printf("\n Insert parameters in microns \n");
3084      printf("\n Insert the x alignment parameter: \n");
3085      cin.getline(value,256);
3086      ixalig = atoi(value);
3087      xalig = (Float_t)ixalig/1000.;
3088      fwrite(&xalig,sizeof(xalig),1,f);
3089      printf(" xalig = %f mm\n",xalig);
3090      //
3091      printf("\n Insert the y alignment parameter: \n");
3092      cin.getline(value,256);
3093      iyalig = atoi(value);
3094      yalig = (Float_t)iyalig/1000.;
3095      fwrite(&yalig,sizeof(yalig),1,f);
3096      printf(" yalig = %f mm\n",yalig);
3097      //
3098      printf("\n Insert the z alignment parameter: \n");
3099      cin.getline(value,256);
3100      izalig = atoi(value);
3101      zalig = (Float_t)izalig/1000.;
3102      fwrite(&zalig,sizeof(zalig),1,f);
3103      printf(" zalig = %f mm\n",zalig);
3104      //
3105      printf("\n Insert the signal/noise threshold (times 1000): \n");
3106      cin.getline(value,256);
3107      iemip = atoi(value);
3108      emip = (Float_t)iemip/1000.;
3109      fwrite(&emip,sizeof(emip),1,f);
3110      printf(" emip = %f mip \n",emip);
3111      //
3112      fclose(f);
3113    }
3114    
3115    void FCaloREADALIG(TString txtple){
3116      FILE *f;
3117      f = fopen(txtple.Data(),"rb");
3118      Float_t mip = 0.;
3119      fread(&mip,sizeof(mip),1,f);
3120      printf(" xalig = %f \n",mip);
3121      mip = 0.;
3122      fread(&mip,sizeof(mip),1,f);
3123      printf(" yalig = %f \n",mip);
3124      mip = 0.;
3125      fread(&mip,sizeof(mip),1,f);
3126      printf(" zalig = %f \n",mip);
3127      mip = 0.;
3128      fread(&mip,sizeof(mip),1,f);
3129      printf(" signal threshold = %f \n",mip);
3130      fclose(f);
3131  }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.23