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

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

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

revision 1.10 by mocchiut, Wed Aug 12 14:57:00 2009 UTC revision 1.12 by mocchiut, Mon Aug 17 15:00:36 2009 UTC
# Line 671  CaloLong::CaloLong(PamLevel2 *l2p){   Line 671  CaloLong::CaloLong(PamLevel2 *l2p){  
671    umax = 100.;    umax = 100.;
672    slmax = "";    slmax = "";
673    sumax = "";    sumax = "";
674      xyaverage = true;
675    //    //
676  };  };
677    
# Line 994  void CaloLong::Fit(){ Line 995  void CaloLong::Fit(){
995    this->Fit(false);    this->Fit(false);
996  };  };
997    
998  Float_t CaloLong::Evaluate(TString s, Float_t tmax){  Float_t CaloLong::Evaluate(TString s, Float_t tmax, Float_t X0pl){
999    /* SAMPLE OUTPUT:    /* SAMPLE OUTPUT:
1000       Enter Infix Expression : A + B + C / (E - F)       Enter Infix Expression : A + B + C / (E - F)
1001       Postfix Expression is : A B + C E F - / +         Postfix Expression is : A B + C E F - / +  
1002    */      */  
1003    if ( !s.Contains("t") ){    if ( !s.Contains("tmax") && !s.Contains("X0pl") ){
1004      printf(" ERROR, the input formula must contain \"t\"\n");      printf(" ERROR, the input formula must contain \"t\"\n");
1005      return 0.;      return 0.;
1006    };    };
# Line 1008  Float_t CaloLong::Evaluate(TString s, Fl Line 1009  Float_t CaloLong::Evaluate(TString s, Fl
1009      return 0.;      return 0.;
1010    };    };
1011    TString g=Form("%f",tmax);    TString g=Form("%f",tmax);
1012      TString h=Form("%f",X0pl);
1013    TString *ts= new TString("");    TString *ts= new TString("");
1014    ts->Prepend(s.Data());    ts->Prepend(s.Data());
1015    ts->ReplaceAll("t",1,g.Data(),g.Capacity());    ts->ReplaceAll("tmax",4,g.Data(),g.Capacity());
1016      ts->ReplaceAll("X0pl",4,h.Data(),h.Capacity());
1017    ts->Prepend("(");    ts->Prepend("(");
1018    ts->Append(")");    ts->Append(")");
1019    if ( debug )  printf(" ts %s tssize %i capac %i s %s g %s \n",ts->Data(),ts->Sizeof(),ts->Capacity(),s.Data(),g.Data());    if ( debug )  printf(" ts %s tssize %i capac %i s %s g %s \n",ts->Data(),ts->Sizeof(),ts->Capacity(),s.Data(),g.Data());
# Line 1074  void CaloLong::Fit(Bool_t draw){ Line 1077  void CaloLong::Fit(Bool_t draw){
1077    };    };
1078    //    //
1079    TString thid = Form("hlongfit");          TString thid = Form("hlongfit");      
1080    TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));    TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
1081    if ( th ) th->Delete();    if ( th ) th->Delete();
1082      //
1083      TString ghid = Form("glongfit");      
1084      TGraphErrors *gh  = dynamic_cast<TGraphErrors*>(gDirectory->FindObject(ghid));
1085      if ( gh ) gh->Delete();
1086      //
1087    Float_t xpos = 0.;    Float_t xpos = 0.;
1088    Float_t enemip = 0.;    Float_t enemip = 0.;
1089    Float_t xmax = NC * X0pl + 0.2;    Float_t xmax = NC * X0pl + 0.2;
1090    //  th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax);    //
1091    th = new TH1F(thid,thid,100,-0.2,xmax);    Double_t xxx[50];
1092      Double_t exx[50];
1093      Double_t yyy[50];
1094      Double_t eyy[50];
1095      Int_t numpo = 0;
1096      memset(xxx,0,50*sizeof(Double_t));
1097      memset(exx,0,50*sizeof(Double_t));
1098      memset(yyy,0,50*sizeof(Double_t));
1099      memset(eyy,0,50*sizeof(Double_t));
1100    //    //
1101    // AGH, BUG!    // AGH, BUG!
1102    //    //
# Line 1094  void CaloLong::Fit(Bool_t draw){ Line 1110  void CaloLong::Fit(Bool_t draw){
1110      mmax = NC;      mmax = NC;
1111    };    };
1112    //    //
1113      Float_t emax = 0.;
1114    Float_t qtotparz = 0.;    Float_t qtotparz = 0.;
1115    for (Int_t st=mmin;st<mmax+1;st++){    for (Int_t st=mmin;st<mmax+1;st++){
1116      enemip = 0.;      enemip = 0.;
1117      xpos = (st - mmin) * X0pl;      xpos = (st - mmin) * X0pl;
1118      if ( st > mmin && st < mmax ){            //
1119        if ( no18x && ( st == 18+1 || st == mask18b+1 )){      if ( xyaverage ){
1120          if ( !maskYO ){        //
1121            enemip = 2. * eplane[1][st];        if ( st > mmin && st < mmax ){      
1122            if ( no18x && ( st == 18+1 || st == mask18b+1 )){
1123              if ( !maskYO ){
1124                enemip = 2. * eplane[1][st];
1125              } else {
1126                enemip = eplane[1][st];
1127              };
1128          } else {          } else {
1129            enemip = eplane[1][st];            enemip = eplane[0][st-1] + eplane[1][st];
1130          };          };
1131        } else {        } else {
1132          enemip = eplane[0][st-1] + eplane[1][st];          if ( st == mmin ){
1133              if ( !maskYE ){
1134                enemip = 2. * eplane[1][st];
1135              } else {
1136                enemip = eplane[1][st];
1137              };
1138            };
1139            if ( st == mmax ){
1140              if ( !maskXE ){
1141                enemip = 2. * eplane[0][st-1];
1142              } else {
1143                enemip = eplane[0][st-1];
1144              };
1145            };
1146          };
1147          //
1148          qtotparz += enemip;
1149          if ( enemip > emax ) emax = enemip;
1150          if ( enemip > 0. ){
1151            xxx[numpo] = xpos;
1152            exx[numpo] = 0.1;
1153            yyy[numpo] = enemip;
1154            eyy[numpo] = sqrt(enemip*3.)+sqrt(5.);
1155            numpo++;
1156            //      th->Fill(xpos,enemip);
1157            if ( debug ) printf(" AVE Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
1158        };        };
1159      } else {      } else {
1160        if ( st == mmin ){        for (Int_t jj=0; jj<2; jj++){
1161          if ( !maskYE ){          if ( st > mmin && st < mmax ){
1162            enemip = 2. * eplane[1][st];            if ( jj == 0 && no18x && ( st == 18+1 || st == mask18b+1 )){
1163                enemip = eplane[1][st];
1164              } else {
1165                if ( jj == 0 ){
1166                  enemip = eplane[jj][st-1];
1167                } else {
1168                  enemip = eplane[jj][st];
1169                };        
1170              };
1171          } else {          } else {
1172            enemip = eplane[1][st];            if ( st == mmin && jj == 1 ){
1173                enemip = eplane[jj][st];
1174              };
1175              if ( st == mmax && jj == 0){
1176                enemip = eplane[jj][st-1];
1177              };
1178          };          };
1179        };          //
1180        if ( st == mmax ){          qtotparz += enemip;
1181          if ( !maskXE ){        if ( enemip > emax ) emax = enemip;
1182            enemip = 2. * eplane[0][st-1];          if ( enemip > 0. ){
1183          } else {            xxx[numpo] = xpos;
1184            enemip = eplane[0][st-1];            exx[numpo] = 0.1;
1185              yyy[numpo] = enemip;
1186              eyy[numpo] = sqrt(enemip*3.)+sqrt(5.);
1187              //      eyy[numpo] = sqrt(enemip)/(st*0.95);
1188              numpo++;
1189              //      th->Fill(xpos,enemip);
1190              if ( debug ) printf(" XY Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
1191          };          };
1192        };        };        
     };  
     //  
     qtotparz += enemip;  
     if ( enemip > 0. ){  
       th->Fill(xpos,enemip);  
       if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip);  
1193      };      };
1194    
1195      //      //
1196      //    for (Int_t v=1; v>=0;v--)// {      //    for (Int_t v=1; v>=0;v--)// {
1197  //       //  //       //
# Line 1155  void CaloLong::Fit(Bool_t draw){ Line 1217  void CaloLong::Fit(Bool_t draw){
1217  //     };  //     };
1218    };    };
1219    //    //
1220    TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);    //  th = new TH2F(thid,thid,int(NC*1.5),-0.2,xmax);
1221      th = new TH2F(thid,thid,1000,-0.2,xmax,1000,0.,emax*1.2);
1222      gh = new TGraphErrors(numpo,xxx,yyy,exx,eyy);
1223      TF1 *lfit  = dynamic_cast<TF1*>(gDirectory->FindObject("lfit"));
1224      if ( lfit ) lfit->Delete();
1225      lfit = new TF1("lfit",ccurve,0.,xmax,3);
1226    if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz);    if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz);
1227    E0 = qtotparz;    E0 = qtotparz;
1228    //  E0 = clp->qtot;    //  E0 = clp->qtot;
# Line 1166  void CaloLong::Fit(Bool_t draw){ Line 1233  void CaloLong::Fit(Bool_t draw){
1233    //  lfit->SetParLimits(0,0.,1000.);    //  lfit->SetParLimits(0,0.,1000.);
1234    //  lfit->SetParLimits(1,-1.,80.);    //  lfit->SetParLimits(1,-1.,80.);
1235    //  lfit->SetParLimits(2,-1.,10.);    //  lfit->SetParLimits(2,-1.,10.);
1236    TString optio;    //  TString optio = "ROW"; // "RO"
1237      TString optio = "RO"; // "RO"
1238    if ( debug ){    if ( debug ){
1239      //    optio = "MERBOV";      optio += "NV";
1240      //    optio = "MEROV";      if ( draw ) optio += "V";
1241      //    optio = "EROV";    } else {
1242      optio = "RNOV";      optio += "NQ";
1243      if ( draw ) optio = "ROV";      if ( draw ) optio += "Q";
   } else {  
     //    optio = "MERNOQ";  
     //    optio = "ERNOQ";  
     optio = "RNOQ";  
     if ( draw ) optio = "ROQ";  
1244    };    };
1245    //    //
1246    if ( debug ) printf(" OK, start the fitting procedure...\n");    if ( debug ) printf(" OK, start the fitting procedure...\n");
1247    //    //
1248    fitresult = th->Fit("lfit",optio);    //  fitresult = th->Fit("lfit",optio);
1249      fitresult = gh->Fit("lfit",optio);
1250    //    //
1251    if ( debug ) printf(" the fit is done! result: %i \n",fitresult);    if ( debug ) printf(" the fit is done! result: %i \n",fitresult);
1252    //    //
# Line 1203  void CaloLong::Fit(Bool_t draw){ Line 1267  void CaloLong::Fit(Bool_t draw){
1267      asymm = -1.;      asymm = -1.;
1268      defE0 = -1.;      defE0 = -1.;
1269    } else {    } else {
1270      if ( slmax.MaybeRegexp() ) lmax = this->Evaluate(slmax,tmax);      if ( slmax.MaybeRegexp() ) lmax = this->Evaluate(slmax,tmax,X0pl);
1271      if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax);      if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax,X0pl);
1272      Int_t npp = 1000;      Int_t npp = 1000;
1273      double *xpp=new double[npp];      double *xpp=new double[npp];
1274      double *wpp=new double[npp];      double *wpp=new double[npp];
# Line 1259  void CaloLong::Fit(Bool_t draw){ Line 1323  void CaloLong::Fit(Bool_t draw){
1323      // axis titles      // axis titles
1324      th->SetXTitle("Depth [X0]");      th->SetXTitle("Depth [X0]");
1325      th->SetYTitle("Energy [MIP]");      th->SetYTitle("Energy [MIP]");
1326      th->DrawCopy("Perror");      //    th->DrawCopy("Perror");
1327        th->DrawCopy();
1328        gh->SetMarkerSize(1.);
1329        gh->SetMarkerStyle(21);
1330        //    gh->SetMarkerColor(kRed);
1331        gh->Draw("Psame");
1332      lfit->Draw("same");      lfit->Draw("same");
1333      tc->Modified();      tc->Modified();
1334      tc->Update();      tc->Update();
# Line 1271  void CaloLong::Fit(Bool_t draw){ Line 1340  void CaloLong::Fit(Bool_t draw){
1340      if ( th ) th->Delete();      if ( th ) th->Delete();
1341    };    };
1342    //    //
1343    delete lfit;    //  delete lfit;
1344    //    //
1345  };  };
1346    
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.23