--- calo/flight/CaloProfile/src/CaloProfile.cpp	2009/08/12 14:43:35	1.9
+++ calo/flight/CaloProfile/src/CaloProfile.cpp	2010/05/13 13:55:37	1.18
@@ -281,7 +281,9 @@
   PKT = 0;
   atime = 0;
   //
+  suf = "";
   debug = false;
+  usepl18x = false;
   //
 };
 
@@ -297,7 +299,9 @@
   PKT = 0;
   atime = 0;
   //
+  suf = "";
   debug = false;
+  usepl18x = false;
   //
 };
 
@@ -375,7 +379,7 @@
   //
   for (Int_t v=minv; v<maxv;v++){
      for (Int_t p=minp; p<maxp;p++){
-	TString hid = Form("clatv%ip%i",v,p);	
+       TString hid = Form("clatv%ip%i%s",v,p,suf.Data());	
         TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
 	if ( tc ){
 //	 tc->Clear();
@@ -383,7 +387,7 @@
 	 tc = new TCanvas(hid,hid);
 	};
 	//
-	TString thid = Form("hlatv%ip%i",v,p);	
+	TString thid = Form("hlatv%ip%i%s",v,p,suf.Data());	
         TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
 	if ( th ) th->Delete();
 //	 th->Clear();
@@ -426,7 +430,7 @@
   gStyle->SetNdivisions(510,"XY");
   //
   for (Int_t p=minp; p<maxp;p++){
-    TString hid = Form("c2dp%i",p);	
+    TString hid = Form("c2dp%i%s",p,suf.Data());	
     TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
     if ( tc ){
       //	 tc->Clear();
@@ -434,7 +438,7 @@
       tc = new TCanvas(hid,hid);
     };
     //
-    TString thid = Form("h2dp%i",p);	
+    TString thid = Form("h2dp%i%s",p,suf.Data());	
     TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
     if ( th ) th->Delete();
     //	 th->Clear();
@@ -525,6 +529,9 @@
   //
   for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
     mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
+    //
+    if ( !usepl18x && view1==0 && plane1==18 ) mip1 = 0.;
+    //
     estrip[view1][plane1][strip1] = mip1;
   };
   //
@@ -572,6 +579,9 @@
   //
   for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
     mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
+    //
+    if ( !usepl18x && view1==0 && plane1==18 ) mip1 = 0.;
+    //
     es[view1][plane1][strip1] = mip1;
   };
   //
@@ -661,6 +671,7 @@
   mask18b = -1;
   //
   no18x = true;
+  usepl18x = !no18x;
   debug = false;
   maskXE = false;
   maskXO = false;
@@ -671,9 +682,23 @@
   umax = 100.;
   slmax = "";
   sumax = "";
+  suf = "";
+  xyaverage = true;
+  //
+  letmax = -1;
+  heavytail = false;
+  //  lmipth = 100.;
+  lmipth = 0.;
   //
 };
 
+TF1 *CaloLong::GetFit(){
+  TString fnam=Form("lfit%s",suf.Data());
+  TF1 *lfit  = dynamic_cast<TF1*>(gDirectory->FindObject(fnam));
+  if ( lfit ) return lfit;
+  return NULL;
+}
+
 void CaloLong::MaskSection(TString sec){
   sec.ToUpper();
   if ( sec.Contains("XO") ) maskXO = true;
@@ -696,6 +721,8 @@
 
 void CaloLong::Clear(){
   //
+  if ( debug ) printf(" Clear called \n");
+  //
   memset(eplane,0, 2*22*sizeof(Float_t));
   //
   chi2 = 0.;
@@ -838,6 +865,9 @@
   Bool_t gof = true;
   for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){
     mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
+    //
+    if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
+    //
     gof = true;
     if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
     if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
@@ -848,8 +878,8 @@
   //
   // inclination factor (stolen from Daniele's code)
   //
-  Float_t ytgx = 0;
-  Float_t ytgy = 0;
+  Float_t ytgx = 0.;
+  Float_t ytgy = 0.;
   ytgx = 0.76 * clp->tanx[0];
   ytgy = 0.76 * clp->tany[0];  
   X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
@@ -943,13 +973,16 @@
     };
   };
   //
+  if ( !clp && debug ) printf(" ERROR!! no CaloLevel2 obj!!\n");
+  //
   // inclination factor (stolen from Daniele's code)
   //
-  Float_t ytgx = 0;
-  Float_t ytgy = 0;
+  Float_t ytgx = 0.;
+  Float_t ytgy = 0.;
   ytgx = 0.76 * clp->tanx[0];
   ytgy = 0.76 * clp->tany[0];  
   X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
+  if ( debug ) printf(" X0pl %f tanx %f tany %f \n",X0pl,ytgx,ytgy);
   //
   // Find experimental plane of maximum
   //
@@ -970,6 +1003,7 @@
   //
   if ( vmax == 0 ) pmax++;
   etmax = pmax * X0pl;
+  letmax = etmax;
   //
   if ( debug ) this->Print();
   if ( debug ) printf(" exit \n");
@@ -994,19 +1028,25 @@
   this->Fit(false);
 };
 
-Float_t CaloLong::Evaluate(TString s, Float_t tmax){
+Float_t CaloLong::Evaluate(TString s, Float_t tmax, Float_t X0pl){
   /* SAMPLE OUTPUT:
      Enter Infix Expression : A + B + C / (E - F)
      Postfix Expression is : A B + C E F - / +  
   */  
-  if ( !s.Contains("t") ){
+  if ( !s.Contains("tmax") && !s.Contains("X0pl") ){
     printf(" ERROR, the input formula must contain \"t\"\n");
     return 0.;
   };
+  if ( tmax != tmax ){
+    printf(" ERROR, tmax is nan! \n");
+    return 0.;
+  };
   TString g=Form("%f",tmax);
+  TString h=Form("%f",X0pl);
   TString *ts= new TString("");
   ts->Prepend(s.Data());
-  ts->ReplaceAll("t",1,g.Data(),g.Capacity());
+  ts->ReplaceAll("tmax",4,g.Data(),g.Capacity());
+  ts->ReplaceAll("X0pl",4,h.Data(),h.Capacity());
   ts->Prepend("(");
   ts->Append(")");
   if ( debug )  printf(" ts %s tssize %i capac %i s %s g %s \n",ts->Data(),ts->Sizeof(),ts->Capacity(),s.Data(),g.Data());
@@ -1059,7 +1099,7 @@
     gStyle->SetNdivisions(510,"XY");
   };
   //
-  TString hid = Form("clongfit");	
+  TString hid = Form("clongfit%s",suf.Data());	
   TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
   //  if ( tc ) tc->Delete();
   //  if ( tc ) tc->Close();
@@ -1069,14 +1109,27 @@
     if ( tc ) tc->cd();
   };
   //
-  TString thid = Form("hlongfit");	
-  TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
+  TString thid = Form("hlongfit%s",suf.Data());	
+  TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
   if ( th ) th->Delete();
+  //
+  TString ghid = Form("glongfit%s",suf.Data());	
+  TGraphErrors *gh  = dynamic_cast<TGraphErrors*>(gDirectory->FindObject(ghid));
+  if ( gh ) gh->Delete();
+  //
   Float_t xpos = 0.;
   Float_t enemip = 0.;
   Float_t xmax = NC * X0pl + 0.2;
-  //  th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax);
-  th = new TH1F(thid,thid,100,-0.2,xmax);
+  //
+  Double_t xxx[50];
+  Double_t exx[50];
+  Double_t yyy[50];
+  Double_t eyy[50];
+  Int_t numpo = 0;
+  memset(xxx,0,50*sizeof(Double_t));
+  memset(exx,0,50*sizeof(Double_t));
+  memset(yyy,0,50*sizeof(Double_t));
+  memset(eyy,0,50*sizeof(Double_t));
   //
   // AGH, BUG!
   //
@@ -1090,42 +1143,99 @@
     mmax = NC;
   };
   //
+  Float_t emax = 0.;
   Float_t qtotparz = 0.;
   for (Int_t st=mmin;st<mmax+1;st++){
     enemip = 0.;
     xpos = (st - mmin) * X0pl;
-    if ( st > mmin && st < mmax ){      
-      if ( no18x && ( st == 18+1 || st == mask18b+1 )){
-	if ( !maskYO ){
-	  enemip = 2. * eplane[1][st]; 
+    //
+    if ( xyaverage ){
+      //
+      if ( st > mmin && st < mmax ){      
+	if ( no18x && ( st == 18+1 || st == mask18b+1 )){
+	  if ( !maskYO ){
+	    enemip = 2. * eplane[1][st]; 
+	  } else {
+	    enemip = eplane[1][st]; 
+	  };
 	} else {
-	  enemip = eplane[1][st]; 
+	  enemip = eplane[0][st-1] + eplane[1][st]; 
 	};
       } else {
-	enemip = eplane[0][st-1] + eplane[1][st]; 
+	if ( st == mmin ){
+	  if ( !maskYE ){
+	    enemip = 2. * eplane[1][st]; 
+	  } else {
+	    enemip = eplane[1][st]; 
+	  };
+	};
+	if ( st == mmax ){
+	  if ( !maskXE ){
+	    enemip = 2. * eplane[0][st-1]; 
+	  } else {
+	    enemip = eplane[0][st-1]; 
+	  };
+	};
+      };
+      //
+      qtotparz += enemip;
+      if ( enemip > emax ) emax = enemip;
+      if ( enemip > 0. ){
+	xxx[numpo] = xpos;
+	exx[numpo] = 0.1;
+	yyy[numpo] = enemip;
+	eyy[numpo] = sqrt(enemip*3.)+sqrt(5.);
+	//if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo;
+	if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = sqrt(enemip)/5.;
+	if ( xpos > letmax-1 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 5.;
+	//if ( xpos > letmax-2 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 7.;
+	if ( xpos < 3. && heavytail ) eyy[numpo] /= 5.;
+	numpo++;
+	//	th->Fill(xpos,enemip);
+	if ( debug ) printf(" AVE Filling: st %i xpos %f energy %f error %f letmax %f lmpith %f \n",st,xpos,enemip,eyy[numpo-1],letmax,lmipth);
       };
     } else {
-      if ( st == mmin ){
-	if ( !maskYE ){
-	  enemip = 2. * eplane[1][st]; 
+      for (Int_t jj=0; jj<2; jj++){
+	if ( st > mmin && st < mmax ){ 
+	  if ( jj == 0 && no18x && ( st == 18+1 || st == mask18b+1 )){
+	    enemip = eplane[1][st]; 
+	  } else {
+	    if ( jj == 0 ){
+	      enemip = eplane[jj][st-1];
+	    } else {
+	      enemip = eplane[jj][st]; 
+	    };	      
+	  };
 	} else {
-	  enemip = eplane[1][st]; 
+	  if ( st == mmin && jj == 1 ){
+	    enemip = eplane[jj][st]; 
+	  };
+	  if ( st == mmax && jj == 0){
+	    enemip = eplane[jj][st-1]; 
+	  };
 	};
-      };
-      if ( st == mmax ){
-	if ( !maskXE ){
-	  enemip = 2. * eplane[0][st-1]; 
-	} else {
-	  enemip = eplane[0][st-1]; 
+	//
+	qtotparz += enemip;
+      if ( enemip > emax ) emax = enemip;
+	if ( enemip > 0. ){
+	  xxx[numpo] = xpos;
+	  exx[numpo] = 0.1;
+	  yyy[numpo] = enemip;
+	  eyy[numpo] = sqrt(enemip*3.)+sqrt(5.);
+	  //if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo;
+	  //if ( xpos > letmax && enemip > lmipth && heavytail ) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo;
+	  if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = sqrt(enemip)/5.;
+	  if ( xpos > letmax-1 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 5.;
+	//if ( xpos > letmax-2 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 7.;
+	  if ( xpos < 3. && heavytail ) eyy[numpo] /= 5.;	 
+	  //	  eyy[numpo] = sqrt(enemip)/(st*0.95);
+	  numpo++;
+	  //	  th->Fill(xpos,enemip);
+	  if ( debug ) printf(" XY Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
 	};
-      };
-    };
-    //
-    qtotparz += enemip;
-    if ( enemip > 0. ){
-      th->Fill(xpos,enemip);
-      if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
+      };	
     };
+
     //
     //    for (Int_t v=1; v>=0;v--)// {
 //       //
@@ -1151,7 +1261,21 @@
 //     };
   };
   //
-  TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);
+  //  th = new TH2F(thid,thid,int(NC*1.5),-0.2,xmax);
+  th = new TH2F(thid,thid,1000,-0.2,xmax,1000,0.,emax*1.2);
+  gh = new TGraphErrors(numpo,xxx,yyy,exx,eyy);
+  TString fnam=Form("lfit%s",suf.Data());
+  TF1 *lfit  = dynamic_cast<TF1*>(gDirectory->FindObject(fnam));
+  if ( lfit ) lfit->Delete();
+  lfit = new TF1(fnam,ccurve,0.,xmax,3);
+//   if ( !lfit ){
+//     lfit = new TF1(fnam,ccurve,0.,xmax,3);
+//   };
+//   lfit->Clear();
+//   lfit->SetName(fnam);
+//   lfit->SetRange(0.,xmax);
+//  lfit->Print();
+  //
   if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz); 
   E0 = qtotparz;
   //  E0 = clp->qtot;
@@ -1162,23 +1286,22 @@
   //  lfit->SetParLimits(0,0.,1000.);
   //  lfit->SetParLimits(1,-1.,80.);
   //  lfit->SetParLimits(2,-1.,10.);
-  TString optio;
+  // TString optio = "ROW"; // "RO"
+  TString optio = "RO"; // "RO"
   if ( debug ){
-    //    optio = "MERBOV";
-    //    optio = "MEROV";
-    //    optio = "EROV";
-    optio = "RNOV";
-    if ( draw ) optio = "ROV";
+    optio += "NV";
+    if ( draw ) optio += "V";
   } else {
-    //    optio = "MERNOQ";
-    //    optio = "ERNOQ";
-    optio = "RNOQ";
-    if ( draw ) optio = "ROQ";
+    optio += "NQ";
+    if ( draw ) optio += "Q";
   };
   //
-  if ( debug ) printf(" OK, start the fitting procedure...\n");
+  if ( debug ) printf(" OK, start the fitting procedure... options: %s \n",optio.Data());
   //
-  fitresult = th->Fit("lfit",optio);
+  //  fitresult = th->Fit("lfit",optio);
+  //  lfit->Update();
+  fitresult = gh->Fit("lfit",optio);
+  //  fitresult = gh->Fit("lfit","ROQW");
   //
   if ( debug ) printf(" the fit is done! result: %i \n",fitresult);
   //
@@ -1199,8 +1322,8 @@
     asymm = -1.;
     defE0 = -1.;
   } else {
-    if ( slmax.MaybeRegexp() ) lmax = this->Evaluate(slmax,tmax);
-    if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax);
+    if ( slmax.MaybeRegexp() ) lmax = this->Evaluate(slmax,tmax,X0pl);
+    if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax,X0pl);
     Int_t npp = 1000;
     double *xpp=new double[npp];
     double *wpp=new double[npp];
@@ -1222,6 +1345,8 @@
     delete w;
     delete xp;
     delete wp;
+    delete xpp;
+    delete wpp;
     //    Float_t i10max = lfit->Integral(0.,10.*tmax);
     if ( debug ) printf(" Integral: %f \n",i10max); 
     //
@@ -1255,7 +1380,12 @@
     // axis titles
     th->SetXTitle("Depth [X0]");
     th->SetYTitle("Energy [MIP]");
-    th->DrawCopy("Perror");
+    //    th->DrawCopy("Perror");
+    th->DrawCopy();
+    gh->SetMarkerSize(1.);
+    gh->SetMarkerStyle(21);
+    //    gh->SetMarkerColor(kRed);
+    gh->Draw("Psame");
     lfit->Draw("same");
     tc->Modified();
     tc->Update();
@@ -1267,7 +1397,16 @@
     if ( th ) th->Delete();
   };
   //
-  delete lfit;
+//   gMinuit->SetPrintLevel(1);
+//   gMinuit->mnamin();
+//   //  gMinuit->mnfree(0);
+//   gMinuit->mncler();
+//   gMinuit->mnrset(1);
+  //  gMinuit->Delete();
+  //ROOT::Math::Minimizer::Clear();
+  //  gMinuit->mnhelp("*");
+  //
+  //  delete lfit;
   //
 };
 
@@ -1296,7 +1435,7 @@
   //
   if ( maxv != -1 ){
     for (Int_t v=minv; v<maxv;v++){
-      TString hid = Form("clongv%i",v);	
+      TString hid = Form("clongv%i%s",v,suf.Data());	
       TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
       if ( tc ){
 	//	 tc->Clear();
@@ -1304,7 +1443,7 @@
 	tc = new TCanvas(hid,hid);
       };
       //
-      TString thid = Form("hlongv%i",v);	
+      TString thid = Form("hlongv%i%s",v,suf.Data());	
       TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
       if ( th ) th->Delete();
       //	 th->Clear();
@@ -1323,14 +1462,14 @@
     };
   } else {
     //
-    TString hid = Form("clongvyvx");	
+    TString hid = Form("clongvyvx%s",suf.Data());	
     TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
     if ( tc ){
     } else {
       tc = new TCanvas(hid,hid);
     };
     //
-    TString thid = Form("hlongvyvx");	
+    TString thid = Form("hlongvyvx%s",suf.Data());	
     TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
     if ( th ) th->Delete();
     th = new TH1F(thid,thid,44,-0.5,43.5);