diff --git a/.gitignore b/.gitignore
index 39cbaf6..c0e4fb9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,6 +3,7 @@
 *.so
 *.root
 *.evt
+*.pdf
 
 root_data
 EventBuilder
diff --git a/DecayFinder.C b/DecayFinder.C
index 139ba6e..203dba2 100644
--- a/DecayFinder.C
+++ b/DecayFinder.C
@@ -13,11 +13,11 @@
 //===================== Settings
 Bool_t debug = false;
 
-TString pidCutFileName = "PIDCuts.root";
+TString pidCutFileName = ""; //"PIDCuts.root";
 
-double distThreshold = 0.2; 
-double forwardTime = 250; //ms
-double backwardTime = 100; //ms
+double distThreshold = 0.1; 
+double forwardTime = 100; //ms
+double backwardTime = 50; //ms
 
 //===================== histograms
 
@@ -31,7 +31,7 @@ TH2F * hBetaPos;
 
 //===================== Parameters
 
-int tick2ns = 8;
+int tick2ns = 4;
 double tick2ms = tick2ns / 1e6;
 
 TFile * cutFile;
@@ -41,12 +41,44 @@ int numCut = 0;
 
 int numMatch = 0;
 
+vector<vector<int>> isotopes = {
+                            {19,  7}, {20,  7}, {21,  7},
+                            {22,  8}, {23,  8}, {24,  8},
+                  {24,  9}, {25,  9}, {26,  9}, {27,  9},           {29,  9},
+                  {27, 10}, {28, 10}, {29, 10}, {30, 10}, {31, 10}, {32, 10},
+        {29, 11}, {30, 11}, {31, 11}, {32, 11}, {33, 11}, {34, 11}, {35, 11},
+        {32, 12}, {33, 12}, {34, 12}, {35, 12}, {36, 12}, {37, 12}, {38, 12},
+        {35, 13}, {36, 13}, {37, 13}, {38, 13}, {39, 13}, {40, 13}, {41, 13},
+                  {39, 14}, {40, 14}, {41, 14}, {42, 14}, 
+                            {43, 15}, {44, 15}
+                          }; 
+
+string isoSym(int Z){
+  switch (Z) {
+    case  3 : return "Li"; break;
+    case  4 : return "Be"; break;
+    case  5 : return "B"; break;
+    case  6 : return "C"; break;
+    case  7 : return "N"; break;
+    case  8 : return "O"; break;
+    case  9 : return "F"; break;
+    case 10 : return "Ne"; break;
+    case 11 : return "Na"; break;
+    case 12 : return "Mg"; break;
+    case 13 : return "Al"; break;
+    case 14 : return "Si"; break;
+    case 15 : return "P"; break;
+    case 16 : return "S"; break;
+  }
+  return "XX";
+}
+
 //===================== Begins
 void DecayFinder::Begin(TTree * /*tree*/){
 
   TString option = GetOption();
 
-  hPID = new TH2F("hPID", "PID; TOF; dE", 400, -260, -220, 400, 0, 6500);
+  hPID = new TH2F("hPID", "PID (beam & Implant & back veto); A/Q; Z", 400, 2.2, 4.0, 400, 3, 16);
   
   hDist = new TH2F("hDist", "dist; dx; dy", 100, -distThreshold, distThreshold, 100, -distThreshold, distThreshold);
   
@@ -54,10 +86,35 @@ void DecayFinder::Begin(TTree * /*tree*/){
   hBetaPos = new TH2F("hBetaPos", "Beta Pos; x; y", 200, 0, 1, 200, 0, 1);
   
   //-------------- GetCut;
-  cutFile = new TFile(pidCutFileName, "UPDATE");
+  cutFile = new TFile(pidCutFileName, "READ");
   bool listExist = cutFile->GetListOfKeys()->Contains("cutList");
   if( !listExist ) {
     cutList = new TObjArray();
+    numCut = (int) isotopes.size();
+    
+    printf("============= no Cut found, creating new default Cuts\n");
+    
+    for( int i = 0 ; i < numCut ; i++){
+      
+      cut = new TCutG();
+      
+      printf("%2d | A: %2d, Z: %d \n", i,  isotopes[i][0], isotopes[i][1]); 
+        
+      TString name; name.Form("cut%02d%s", isotopes[i][0], isoSym(isotopes[i][1]).c_str());
+      cut->SetName(name);
+      cut->SetVarX("AoQ");
+      cut->SetVarY("Z");
+      cut->SetTitle(Form("cut%2d", i));
+      cut->SetLineColor(i+1);
+      
+      cut->SetPoint(0, (isotopes[i][0] - 0.45)/isotopes[i][1], isotopes[i][1] - 0.45);
+      cut->SetPoint(1, (isotopes[i][0] - 0.45)/isotopes[i][1], isotopes[i][1] + 0.45);
+      cut->SetPoint(2, (isotopes[i][0] + 0.45)/isotopes[i][1], isotopes[i][1] + 0.45);
+      cut->SetPoint(3, (isotopes[i][0] + 0.45)/isotopes[i][1], isotopes[i][1] - 0.45);
+      cut->SetPoint(4, (isotopes[i][0] - 0.45)/isotopes[i][1], isotopes[i][1] - 0.45);
+
+      cutList->AddAtAndExpand(cut, i);
+    }
   }else{
     cutList = (TObjArray*) cutFile->FindObjectAny("cutList");
     numCut = cutList->GetLast()+1;
@@ -75,8 +132,8 @@ void DecayFinder::Begin(TTree * /*tree*/){
   hvetoF = new TH1F *[numCut];
   
   for( int i = 0; i < numCut; i++){
-    hDecay[i] = new TH1F(Form("hDecay%02d",i), Form("Decay cut-%02d ; [ms]; count", i), 50, -backwardTime, forwardTime); 
-    hvetoF[i] = new TH1F(Form("hvetoF%02d",i), Form("veto-F cut-%02d ; [ms]; count", i), 100, 0, 6000); 
+    hDecay[i] = new TH1F(Form("hDecay%02d%s",isotopes[i][0], isoSym(isotopes[i][1]).c_str()), Form("Decay cut-%02d ; [ms]; count", i), 200, -backwardTime, forwardTime); 
+    hvetoF[i] = new TH1F(Form("hvetoF%02d%s",isotopes[i][0], isoSym(isotopes[i][1]).c_str()), Form("veto-F cut-%02d ; [ms]; count", i), 100, 0, 6000); 
   }
   
   printf("=====================================\n");
@@ -90,10 +147,10 @@ double dist(double x1, double x2, double y1, double y2){
 //===================== Process
 Bool_t DecayFinder::Process(Long64_t entry){
   
-  if( entry > (Long64_t)1e6 ) return kTRUE;
+  //if( entry > (Long64_t)1000 ) return kTRUE;
 
-  b_TOF->GetEntry(entry);
-  b_energy->GetEntry(entry);
+  b_Z->GetEntry(entry);
+  b_AoQ->GetEntry(entry);
   b_crossTime->GetEntry(entry);
   b_crossEnergy->GetEntry(entry);
   b_flag->GetEntry(entry);
@@ -111,28 +168,34 @@ Bool_t DecayFinder::Process(Long64_t entry){
   
   if( debug && entry % 10 == 0 ) printf("------------- %llu\n", entry);
   
-  if( flag & 1 ) hPID->Fill(TOF, energy);
+  if( (flag == 3) && (vetoFlag & 2) == 0 ) hPID->Fill(AoQ, Z);
   
   int cutID = -1;
   for(int i = 0; i < numCut; i++){
     cut = (TCutG*) cutList->At(i);
-    if( cut->IsInside(TOF, energy) ) cutID = i;
+    if( cut->IsInside(AoQ, Z) ) cutID = i;
   }
-  
   if( cutID == -1 ) return kTRUE;
-  if( (flag & 1) == 0  ) return kTRUE; /// no beam
-  //if( (flag & 2) == 0  ) return kTRUE; /// no Ions
-  if( (flag & 4) == 0  ) return kTRUE; /// no beta
+
+  if( flag >= 4 ) return kTRUE;
+  ///if( flag != 3 ) return kTRUE; ///
+
+  //if( (flag & 1) == 0 ) return kTRUE; /// no beam
+  //if( (flag & 2) == 0 ) return kTRUE; /// no Ions
+  //if( (flag & 4) == 4 ) return kTRUE; /// has beta
 
   //if( veto_f > 0 ) return kTRUE;
-  if( veto_r > 0 ) return kTRUE;
-  if( dyIonsTime[0] == 0 ) return kTRUE;
+  //if( veto_r > 0 ) return kTRUE;
+  //if( dyIonsTime[0] == 0 ) return kTRUE;
 
   ImplantEntry = entry;
   ImplantTime = dyIonsTime[0];
   ImplantX = xIons;
   ImplantY = yIons;
   
+  if( TMath::IsNaN(xIons) ) return kTRUE;
+  if( TMath::IsNaN(yIons) ) return kTRUE;
+  
   hIonsPos->Fill( xIons, yIons);
     
   if( debug ) printf("===========%8lld, %13llu, %f, %f\n", ImplantEntry, ImplantTime, ImplantX, ImplantY); 
@@ -151,12 +214,13 @@ Bool_t DecayFinder::Process(Long64_t entry){
     if ( k == ImplantEntry ) continue;
     
     //if( crossTime > 0 ) continue;
-    if( crossEnergy > 0 ) continue;
-    //if( (flag & 2) == 2 ) continue; /// has Ions
+    //if( crossEnergy > 0 ) continue;
+    if( (flag & 2) == 2 ) continue; /// has Ions
+    //if( flag  != 4 ) continue; /// no Beta
     if( (flag & 4) == 0 ) continue; /// no Beta
     if( !TMath::IsNaN(veto_f) ) continue; 
     if( !TMath::IsNaN(veto_r) ) continue; 
-    if( dyBetaTime[0] == 0) continue;
+    //if( dyBetaTime[0] == 0) continue;
     if( dist(ImplantX, xBeta, ImplantY, yBeta) > distThreshold ) continue;
     
     if( k < ImplantEntry && ( ImplantTime - dyBetaTime[0] ) * tick2ms > backwardTime) continue;
@@ -165,7 +229,7 @@ Bool_t DecayFinder::Process(Long64_t entry){
     clock.Stop("timer");
     Double_t time = clock.GetRealTime("timer");
     clock.Start("timer");
-    printf(" %llu[%.2f%%], searched next %lld events | match case %d | expect : %5.2f min\r", 
+    printf(" %llu[%4.2f%%], searched next %lld events | match case %d | expect : %5.2f min\r", 
          ImplantEntry, ImplantEntry*100./totNumEntry, k - ImplantEntry, 
          numMatch, totNumEntry*time/(ImplantEntry+1.)/60.);
       
@@ -184,7 +248,7 @@ Bool_t DecayFinder::Process(Long64_t entry){
     }
   
     hDecay[cutID]->Fill(dT);
-    hvetoF[cutID]->Fill(crossEnergy);
+    //hvetoF[cutID]->Fill(crossEnergy);
     hBetaPos->Fill(xBeta, yBeta);
     
   }
@@ -212,12 +276,12 @@ void DecayFinder::Terminate(){
   int padID = 1; cDecay->cd(padID); 
   hPID->Draw("colz");
   for( int i = 0; i < numCut; i++){
-  	cut = (TCutG*) cutList->At(i);
-  	cut->SetLineColor(i+2);
-  	cut->Draw("same");
+    cut = (TCutG*) cutList->At(i);
+    cut->SetLineColor(i+2);
+    cut->Draw("same");
   }
   
-  padID++; cDecay->cd(padID); //cDecay->cd(padID)->SetLogy();
+  padID++; cDecay->cd(padID); cDecay->cd(padID)->SetLogy();
   double yMax = 0;
   for( int i = 0; i < numCut; i++){
     if( hDecay[i]->GetMaximum() > yMax ) yMax = hDecay[i]->GetMaximum();
@@ -225,7 +289,7 @@ void DecayFinder::Terminate(){
   for( int i = 0; i < numCut; i++){
     hDecay[i]->SetLineColor(i+2);
     hDecay[i]->SetMaximum(yMax *1.2);
-  	hDecay[i]->Draw(i == 0 ? "" : "same");
+    hDecay[i]->Draw(i == 0 ? "" : "same");
   }
   
   //padID++; cDecay->cd(padID);
@@ -242,5 +306,20 @@ void DecayFinder::Terminate(){
 
   padID++; cDecay->cd(padID);
   hBetaPos->Draw("colz");
+  
+  //================ Save historgram 
+  TString fHistRootName = "zzz_all_decay.root";
+  if( true ){
+    TFile * fHist = new TFile(fHistRootName, "recreate");
+    fHist->cd();
+    
+    for( int i = 0 ; i < numCut; i++){
+      hDecay[i]->Write();
+    }
+    
+    fHist->Close();
+    
+    printf("---- Save PID histogram as %s \n", fHistRootName.Data());
+  }
 
 }
diff --git a/DecayFinder.h b/DecayFinder.h
index 93aee38..be9e11f 100644
--- a/DecayFinder.h
+++ b/DecayFinder.h
@@ -24,9 +24,11 @@ public :
 
    // Declaration of leaf types
    ULong64_t       eventID;
-   //Int_t           runID;
-   Double_t        TOF;
-   Double_t        energy;
+   Int_t           runID;
+   //Double_t        TOF;
+   //Double_t        energy;
+   Double_t        Z;
+   Double_t        AoQ;
    UInt_t          crossEnergy;
    ULong64_t       crossTime;
    Short_t         flag;
@@ -42,9 +44,11 @@ public :
 
    // List of branches
    TBranch        *b_eventID;   //!
-   //TBranch        *b_runID;   //!
-   TBranch        *b_TOF;   //!
-   TBranch        *b_energy;   //!
+   TBranch        *b_runID;   //!
+   //TBranch        *b_TOF;   //!
+   //TBranch        *b_energy;   //!
+   TBranch        *b_Z;   //!
+   TBranch        *b_AoQ;   //!
    TBranch        *b_crossTime;   //!
    TBranch        *b_crossEnergy;   //!
    TBranch        *b_flag;    //!
@@ -94,9 +98,12 @@ void DecayFinder::Init(TTree *tree){
    fChain->SetMakeClass(1);
 
    fChain->SetBranchAddress("eventID", &eventID, &b_eventID);
-   //fChain->SetBranchAddress("runID", &runID, &b_runID);
-   fChain->SetBranchAddress("TOF", &TOF, &b_TOF);
-   fChain->SetBranchAddress("energy", &energy, &b_energy);
+   fChain->SetBranchAddress("runID", &runID, &b_runID);
+   //fChain->SetBranchAddress("TOF", &TOF, &b_TOF);
+   //fChain->SetBranchAddress("energy", &energy, &b_energy);
+   fChain->SetBranchAddress("Z", &Z, &b_Z);
+   fChain->SetBranchAddress("AoQ", &AoQ, &b_AoQ);
+   fChain->SetBranchAddress("crossTime", &crossTime, &b_crossTime);
    fChain->SetBranchAddress("crossTime", &crossTime, &b_crossTime);
    fChain->SetBranchAddress("crossEnergy", &crossEnergy, &b_crossEnergy);
    fChain->SetBranchAddress("flag",   &flag, &b_flag);
diff --git a/Hist_Style/PID_paper.C b/Hist_Style/PID_paper.C
new file mode 100644
index 0000000..13e5947
--- /dev/null
+++ b/Hist_Style/PID_paper.C
@@ -0,0 +1,71 @@
+void PID_paper(double gamma, int start= 1){
+  
+  TFile * f1 = new TFile("pid_total.root");
+  TH2F * hPID = (TH2F*) f1->Get("hPID2");
+  
+  gStyle->SetOptStat("");  /// no state
+  gStyle->SetOptTitle(0);  /// no title
+  
+  TCanvas * ccc = new TCanvas("ccc", "ccc", 1000, 1000);
+  ///ccc->Divide(2,1);
+  
+  ///========= set Z-color
+  
+  ///---- custom color
+  ///Int_t colors[] = {0, 1, 2, 3, 4, 5, 6}; 
+  ///gStyle->SetPalette((sizeof(colors)/sizeof(Int_t)), colors);
+  
+  
+  ///----- build-in color
+  int color = kBeach;  
+  //int color = kGreenBrownTerrain;  
+  //int color = kColorPrintableOnGrey;
+  //int color = kSunset;
+  
+  ///gStyle->SetPalette(kBlueRedYellow);
+  ///gStyle->SetPalette(kRainBow);
+  gStyle->SetPalette(color);TColor::InvertPalette();
+  
+  ///======== set contours
+    
+  TGraph * ggg = new TGraph(); /// this is for displaying the gamma correction
+  
+  const int nCon = 200;  /// number of contours
+  double contours[nCon];
+  
+  double zMax = hPID->GetMaximum();
+  
+  for( int i = 0; i < nCon; i++){
+    double xpos = i*1.0/nCon;
+    double ypos = TMath::Power(xpos, gamma) *0.95;
+    
+    contours[i] = start + ypos*(zMax - start);
+    ggg->AddPoint(xpos, ypos); 
+  }
+    
+  ccc->cd();ccc->cd()->SetGrid(0,0);
+  hPID->GetYaxis()->SetRangeUser(1.5, 16);
+  hPID->SetContour(nCon, contours);
+  hPID->Draw("col");
+  
+  TEllipse * si42 = new TEllipse(3.01, 14, 0.05, 0.55);
+  si42->SetFillStyle(0);
+  si42->SetLineColor(2);
+  si42->SetLineWidth(3);
+  
+  si42->Draw("same");
+  
+  TLatex text;
+  text.SetTextColor(2);
+  text.SetTextSize(0.05);
+  text.DrawLatex(3.05, 14.2, "^{42}Si");
+
+  ///------ in case want to see the gamma correction
+  ///ccc->cd(2); 
+  ///ggg->Draw("AP*");
+  
+  ///====== save as pdf
+  ccc->SaveAs(Form("haha_%3.1f_%d_%d.pdf", gamma, start, color));
+
+
+}
diff --git a/Hist_Style/ipid.C b/Hist_Style/ipid.C
new file mode 100644
index 0000000..e36b499
--- /dev/null
+++ b/Hist_Style/ipid.C
@@ -0,0 +1,109 @@
+{
+
+
+  TFile * f1 = new TFile("PIDFile.root");
+  TH2F * hPID = (TH2F*) f1->Get("ipid");
+
+  gStyle->SetOptStat("");  /// no state
+  gStyle->SetOptTitle(0);  /// no title
+
+  TCanvas * ccc = new TCanvas("ccc", "ccc", 1000, 1000);
+
+  const int nBinX = hPID->GetNbinsX();
+  const int nBinY = hPID->GetNbinsY();
+  const double xBinWidth = hPID->GetXaxis()->GetBinWidth(1);
+  const double yBinWidth = hPID->GetYaxis()->GetBinWidth(1);
+
+  
+  double skewPar[6] = {10.9179, 
+                       0.00416034,
+                       -233.306,
+                       3607.16,
+                       1538.37,
+                       0.000430609};
+                       
+  const double c = 299.792458; /// mm/ns
+  const double FL = 56423 ;  /// mm 
+  const double Brho = 9.04; /// T.mm   
+  const double ma = 931.5;
+  const double beta0 = Brho * c / TMath::Sqrt( 9*ma*ma + Brho*Brho*c*c); // for AoQ = 3
+  const double tofOffset = 504.69; /// ns
+  
+  
+  TH2F * hPID2 = new TH2F("hPID2", "PID; ToF; dE", 500, -280, -210, 500, 0, 5500); 
+  TH2F * hZAoQ = new TH2F("hZAoQ", "PID; A/Q; Z", 400, 2.0, 4.0, 400, 1, 17); 
+  TH1F * hZ    = new TH1F("hZ", "Z; Z; count", 500, 1, 17); 
+
+
+  for( int i = 1; i < nBinX ; i++){
+    for( int j = 1; j < nBinY ; j++){
+      double xpos = hPID->GetXaxis()->GetBinCenter(i) + (gRandom->Rndm()-0.5) * xBinWidth;
+      double ypos = hPID->GetYaxis()->GetBinCenter(j) + (gRandom->Rndm()-0.5) * yBinWidth;
+      
+      double zpos = hPID->GetBinContent(i,j);
+      
+      
+      ///===== Skew correction
+      double ToF = xpos - (skewPar[0] * exp( - skewPar[1] * ypos ) + skewPar[2] - skewPar[5] * (skewPar[3] - ypos) * exp( - pow( (skewPar[3]-ypos)/skewPar[4],2) ) ) - 234 + 12; 
+      
+      
+      hPID2->Fill( ToF, ypos, zpos);
+      
+      double beta = FL/c/(ToF + tofOffset);
+      double gamma = 1./TMath::Sqrt(1.-beta*beta);
+      
+      double Z = sqrt((ypos + 11.9473) / 23.097) * TMath::Power(beta / beta0, 1.3)  ;
+      double AoQ = c * Brho / gamma/ beta / ma;
+
+      Z = -0.732942 + 1.47488 * Z  - 0.100088 * Z*Z  + 0.009867 * Z*Z*Z - 0.000476 * Z*Z*Z*Z + 8.90519e-6 * Z*Z*Z*Z*Z;
+    
+      hZAoQ->Fill(AoQ, Z, zpos);
+      hZ->Fill(Z, zpos);
+      
+    }
+  }
+
+  
+  double gamma = 4.0;
+  int start = 0;
+  int color = kColorPrintableOnGrey;  
+  gStyle->SetPalette(color);TColor::InvertPalette();
+  
+  
+  double zMax = hZAoQ->GetMaximum();
+  
+  const int nCon = 200;  /// number of contours
+  double contours[nCon];
+  
+  for( int i = 0; i < nCon; i++){
+    double xpos = i*1.0/nCon;
+    double ypos = TMath::Power(xpos, gamma) *0.95;
+    contours[i] = start + ypos*(zMax - start);
+  }
+
+  ///ccc->Divide(2,1);
+  
+  ccc->cd(1); ///ccc->cd(1)->SetLogz();
+  hZAoQ->SetContour(nCon, contours);
+  hZAoQ->Draw("col");
+  
+  TEllipse * si42 = new TEllipse(3.01, 14, 0.05, 0.55);
+  si42->SetFillStyle(0);
+  si42->SetLineColor(2);
+  si42->SetLineWidth(3);
+  
+  si42->Draw("same");
+  
+  TLatex text;
+  text.SetTextColor(2);
+  text.SetTextSize(0.05);
+  text.DrawLatex(3.05, 14.2, "^{42}Si");
+
+  //ccc->cd(2);
+  //hZ->Draw("hist");
+  
+  ccc->SaveAs(Form("iPID_%3.1f_%d_%d.pdf", gamma, start, color));
+
+
+
+}
diff --git a/PIDCutCreator.C b/PIDCutCreator.C
index bc29cdf..065e8c4 100644
--- a/PIDCutCreator.C
+++ b/PIDCutCreator.C
@@ -45,7 +45,7 @@ void PIDCutCreator(TH2F* hist,
       for( int k = 0; k < numCut; k++){
          if( cutList->At(k) != NULL ){
             printf("found a cut at %2d \n", k);
-	    cut = (TCutG*) cutList->At(k);
+            cut = (TCutG*) cutList->At(k);
             cut->Draw("same");
          }else{
             printf("     No cut at %2d \n", k);
@@ -68,10 +68,10 @@ void PIDCutCreator(TH2F* hist,
       break;
     }
 
-    TString name; name.Form("cut%dd", countCut);
+    TString name; name.Form("cut%2d", countCut);
     cut->SetName(name);
-    cut->SetVarX("TOF");
-    cut->SetVarY("dE");
+    cut->SetVarX("AoQ");
+    cut->SetVarY("Z");
     cut->SetTitle(Form("cut%2d", countCut));
     cut->SetLineColor(countCut+1);
     printf(" cut-%d \n", countCut);
diff --git a/armory/AutoFit.C b/armory/AutoFit.C
index 70cb7fe..7d9a38b 100644
--- a/armory/AutoFit.C
+++ b/armory/AutoFit.C
@@ -1897,40 +1897,91 @@ void fitNGaussPolSub(TH1F * hist, int degPol,  TString fitFile = "AutoFit_para.t
 }
 
 //########################################
-//###### fit A * Exp(-t/T) + B
+//###### fit B + A * Exp(-t/Th1) + A * Th2/(Th1-Th2)*( Exp(-t/Th1) - Exp(-t/Th2) )
 //########################################
-void fitDecay(TH1F * hist){
+Double_t decayFunc(Double_t *x, Double_t *par){
+
+  ///"[0]*TMath::Exp(-x/[1])+ [2]* [3]/([1]-[3]) * ( Exp[-x/[1]] - Exp[-x/[3]]) + [4]"  
+  
+  double ln2 = par[5];
+  if ( x[0] < 0 ) {
+    return par[4];
+  }else{
+    double d1 = par[0]*TMath::Exp(-x[0]/par[1]*ln2);
+    double d2 = par[2]*(TMath::Exp(-x[0]/par[1]*ln2) - TMath::Exp(-x[0]/par[3]*ln2)) * par[3] / (par[1]-par[3]);
+    return d1 + d2 +par[4]; 
+  }
+  
+}
+
+void fitDecay(TH1F * hist, double decayConst1 = 0, double decayConst2 = 0, bool halfLive = false){
+  
+  printf("===============================================================\n");
+  printf("================= fit decay + conseqent decay  ================\n");
+  printf(" A1 Exp(-t/T1) + A2 * T2/(T1-T2) * ( Exp(-t/T1) - Exp(-t/T2) ) + B ");
+  printf("   decayConst1 = T1  \n");
+  printf("   decayConst2 = T2 : when 0, auto guess, when < 0, set A2 = 0\n");
+  printf("   halfLive = false by default, can switch on to change T1, T2 to be half-live\n");
+  printf(" \n");
+  printf("================================================================\n");
 
   gStyle->SetOptStat("");
   TCanvas * cFitDecay = NULL;
   if( gROOT->FindObjectAny("cFitDecay") == NULL ){
-    cFitDecay = new TCanvas("cFitDecay", Form("fit A * Exp(-t/T) + B | fitDecay"), 800, 1000);
+    cFitDecay = new TCanvas("cFitDecay", Form("fit consqeunce decay | fitDecay"), 800, 1000);
   }else{
     delete gROOT->FindObjectAny("cFitDecay") ;
-    cFitDecay = new TCanvas("cFitDecay", Form("fit A * Exp(-t/T) + B | fitDecay"), 800, 1000);
+    cFitDecay = new TCanvas("cFitDecay", Form("fit consqeunce decay | fitDecay"), 800, 1000);
   }
   
   cFitDecay->Divide(1, 2);
+
+  double xMin = hist->GetXaxis()->GetXmin();
+  double xMax = hist->GetXaxis()->GetXmax();
   
-  double yMax = hist->GetBinContent(hist->GetMaximumBin());
+  int xBinYMax = hist->GetMaximumBin();
+  double yMax = hist->GetBinContent(xBinYMax);
+  double bg = hist->GetBinContent(1);
   hist->SetMaximum(yMax*1.2);
   cFitDecay->cd(1); hist->Draw(); 
-  double estHalf = hist->GetBinCenter(hist->FindLastBinAbove(yMax/2.)) *1.3;
+  double estHalf1 = (decayConst1 == 0 ? hist->GetBinCenter(hist->FindLastBinAbove((yMax-bg)/2.)) *1.3 : decayConst1);
+  double estHalf2 = (decayConst2 == 0 ? estHalf1 * 3 : decayConst2);
+  
+  if( decayConst2 >= 0 ){
+    printf("estimator : A1 = A2 = %.2f, T1 = %.2f, T2 = %.2f, bg = %.2f \n", yMax, estHalf1, estHalf2, bg);
+  }else{
+    printf("estimator : A1 = %.2f, T1 = %.2f, A2 = 0, T2 = 1, bg = %.2f \n", yMax, estHalf1, bg);
+  }
 
-  //printf("estimator : A = %.2f, T = %.2f \n", yMax, estHalf);
-
-  TF1 * fit = new TF1("fit", "[0]*TMath::Exp(-x/[1]/TMath::Log(2))+[2]");
+  TF1 * fit = new TF1("fit", decayFunc, xMin, xMax, 6 );
   fit->SetLineWidth(3);
   fit->SetLineColor(1);
   fit->SetNpx(1000);
   
   fit->SetParameter(0, yMax );
-  fit->SetParameter(1, estHalf );
-  fit->SetParameter(2,  1);
-   
-  fit->SetParLimits(0,       0, 1e9); 
+  fit->SetParameter(1, estHalf1 );
+  fit->SetParameter(2, yMax );
+  fit->SetParameter(3, estHalf2 );
+  fit->SetParameter(4,  bg);
+  
+  if( halfLive == false ) {
+    fit->FixParameter(5, 1 );
+  }else{
+    fit->FixParameter(5, TMath::Log(2) );
+  }
+  
+  fit->SetParLimits(0,       1, 1e9); 
   fit->SetParLimits(1,       0, 1e9); 
-  fit->SetParLimits(2,       0, 1e9); 
+  fit->SetParLimits(4,       0, 1e9); 
+  
+  ///======= put A2 = 0 for decayConst < 0 
+  if( decayConst2 < 0 ) {
+    fit->FixParameter(2, 0); 
+    fit->FixParameter(3, 1); 
+  }else{
+    fit->SetParLimits(2, 0, 1e9); 
+    fit->SetParLimits(3,       0, 1e9); 
+  }
   
   gStyle->SetOptFit(000000);
   hist->Fit("fit", "q");
@@ -1943,32 +1994,53 @@ void fitDecay(TH1F * hist){
   text.SetTextFont(82);
   text.SetTextSize(0.04);
   text.SetTextColor(1);
-  
-  text.DrawLatex(0.6, 0.8, "A e^{#left(-#frac{t}{T Ln(2)}#right)}+B");
-  text.DrawLatex(0.6, 0.75, Form("A : %.f(%.f)", paraA[0], paraE[0])); 
-  text.DrawLatex(0.6, 0.7, Form("T : %.f(%.f)", paraA[1], paraE[1])); 
-  text.DrawLatex(0.6, 0.65, Form("B : %.f(%.f)", paraA[2], paraE[2])); 
+  if( decayConst2 >= 0 ){
+    text.DrawLatex(0.3, 0.8, "A1 e^{#left(-#frac{t}{T1}#right)} + A2 #frac{T2}{T1-T2} (e^{#left(-#frac{t}{T1}#right)}-e^{#left(-#frac{t}{T2}#right)})+B");
+  }else{
+    text.DrawLatex(0.3, 0.8, "A1 e^{#left(-#frac{t}{T1}#right)} +B");
+  }
+  text.DrawLatex(0.3, 0.75, Form("A1 : %.2f(%.2f)", paraA[0], paraE[0])); 
+  text.DrawLatex(0.6, 0.75, Form("T1 : %.2f(%.2f)", paraA[1], paraE[1])); 
+  if( decayConst2 >= 0 ){
+    text.DrawLatex(0.3, 0.70, Form("A2 : %.2f(%.2f)", paraA[2], paraE[2])); 
+    text.DrawLatex(0.6, 0.70, Form("T2 : %.2f(%.2f)", paraA[3], paraE[3])); 
+  }
+  text.DrawLatex(0.6, 0.65, Form(" B : %.2f(%.2f)", paraA[4], paraE[4])); 
   
   double chi2 = fit->GetChisquare();
   int ndf = fit->GetNDF();
-  text.DrawLatex(0.6, 0.6, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
+  double sigma = TMath::Sqrt(paraA[4]);
+  text.DrawLatex(0.6, 0.60, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf/sigma));
 
   
-  TF1 * fite = new TF1("fite", "[0]*TMath::Exp(-x/[1])" );
+  TF1 * fite = new TF1("fite", "[0]*TMath::Exp(-x/[1])", 0, xMax );
+  fite->SetLineColor(4);
+  fite->SetLineWidth(2);
   fite->SetParameter(0, paraA[0]);
   fite->SetParameter(1, paraA[1]);
   fite->Draw("same");
+
+  if( decayConst2 >= 0 ){
+    TF1 * fitd = new TF1("fitd", "[0]*(TMath::Exp(-x/[1])-TMath::Exp(-x/[2]))*[2]/([1]-[2])", 0, xMax );
+    fitd->SetLineColor(6);
+    fitd->SetLineWidth(2);
+    fitd->SetParameter(0, paraA[2]);
+    fitd->SetParameter(1, paraA[1]);
+    fitd->SetParameter(2, paraA[3]);
+    fitd->Draw("same");
+  }
   
-  TF1 * fitbg = new TF1("fitbg", "pol0" );
-  fitbg->SetParameter(0, paraA[2]);
+  TF1 * fitbg = new TF1("fitbg", "pol0", xMin, xMax );
+  fitbg->SetLineColor(1);
+  fitbg->SetLineWidth(1);
+  fitbg->SetParameter(0, paraA[4]);
   fitbg->Draw("same");
   
   cFitDecay->cd(2); 
   PlotResidual(hist, fit);
-
-
 }
 
+
 //#################################################
 //#### fit N-Gauss with pol-n BG using mouse click
 //#################################################
diff --git a/peachCake.C b/peachCake.C
index f1acdb9..6a2ef5b 100644
--- a/peachCake.C
+++ b/peachCake.C
@@ -40,10 +40,10 @@ TH1F * hDecay_veto;
 TH1F * hFlag;
 TH1F * hvetoFlag;
 
-TH1F * hZ3;
-TH1F * hZ6;
-TH1F * hZ10;
-TH1F * hZ14;
+///TH1F * hZ3;
+///TH1F * hZ6;
+///TH1F * hZ10;
+///TH1F * hZ14;
 
 //============ parameters
 //TString cutFileName="PIDCuts.root";
@@ -54,9 +54,8 @@ TH1F * hZ14;
 
 vector<vector<double>> eCorr;
 
-double TOFCorrection(double x){
-   
-  // for run 250
+double TOFCorrection(double x){   
+  ///========= for run 250
   double par[6] = {10.9179, 
                    0.00416034,
                    -233.306,
@@ -64,7 +63,6 @@ double TOFCorrection(double x){
                    1538.37,
                    0.000430609};
   return par[0] * exp( - par[1] * x ) + par[2] - par[5] * (par[3] - x) * exp( - pow( (par[3]-x)/par[4],2) ); 
-  
 }
 
 
@@ -83,15 +81,14 @@ void peachCake::Begin(TTree * /*tree*/){
 
   hZ = createTH1F("hZ", "Z; Z; count", 500, 0, 16);
 
-  hImplantIons = createTH2F("hImplantIons", "Implat low gain; x ; y", 400, 0, 1, 400, 0, 1);
-  hImplantBeta = createTH2F("hImplantBeta", "Implat high gain; x ; y", 400, 0, 1, 400, 0, 1);   
+  hImplantIons = createTH2F("hImplantIons", "Implant low gain; x ; y", 400, 0, 1, 400, 0, 1);
+  hImplantBeta = createTH2F("hImplantBeta", "Implant high gain; x ; y", 400, 0, 1, 400, 0, 1);   
 
-  hImplantX = createTH2F("hImplantX", "Implat X; low ; high", 400, 0, 1, 400, 0, 1);
-  hImplantY = createTH2F("hImplantY", "Implat Y; low ; high", 400, 0, 1, 400, 0, 1);   
+  hImplantX = createTH2F("hImplantX", "Implant X; low ; high", 400, 0, 1, 400, 0, 1);
+  hImplantY = createTH2F("hImplantY", "Implant Y; low ; high", 400, 0, 1, 400, 0, 1);   
 
-
-  hImplantDyIons = createTH2F("hImplantDyIonsow", "Implat low; sum_a; dy", 400, 0, 80000, 400, 0, 80000);
-  hImplantDyBeta = createTH2F("hImplantDyBetaigh", "Implat high; sum_a; dy", 400, 0, 8000, 400, 0, 8000);
+  hImplantDyIons = createTH2F("hImplantDyIonslow", "Implant low; sum_a; dy", 400, 0, 16000, 400, 0, 16000);
+  hImplantDyBeta = createTH2F("hImplantDyBetahigh", "Implant high; sum_a; dy", 400, 0, 80000, 400, 0, 80000);
 
   hDecay = createTH1F("hDecay", "Decay (beta.time - ion.time) ; [ticks]; counts", 100, 0, 2000);
   hDecay_veto = createTH1F("hDecay_veto", "Decay (beta.time - ion.time) [veto]; [ticks]; counts", 100, 0, 2000);
@@ -99,10 +96,10 @@ void peachCake::Begin(TTree * /*tree*/){
   hFlag = createTH1F("hFlag", "Flag ( 1 = beam, 2 = Ions, 4 = Beta)", 8, 0, 8);
   hvetoFlag = createTH1F("hvetoFlag", "vetoFlag ( 1 = front, 4 = rear)", 8, 0, 8);
   
-  hZ3  = createTH1F("hZ3", "Z=3; A ; count ", 200, 6, 12);
-  hZ6  = createTH1F("hZ6", "Z=6; A ; count ", 200, 14, 21);
-  hZ10 = createTH1F("hZ10", "Z=10; A ; count ", 200, 25, 33);
-  hZ14 = createTH1F("hZ14", "Z=14; A ; count ", 200, 38, 44);
+  ///hZ3  = createTH1F("hZ3", "Z=3; A ; count ", 200, 6, 12);
+  ///hZ6  = createTH1F("hZ6", "Z=6; A ; count ", 200, 14, 21);
+  ///hZ10 = createTH1F("hZ10", "Z=10; A ; count ", 200, 25, 33);
+  ///hZ14 = createTH1F("hZ14", "Z=14; A ; count ", 200, 38, 44);
 
   eCorr = LoadCorrectionParameters("correction_e.dat");
   
@@ -127,14 +124,18 @@ Bool_t peachCake::Process(Long64_t entry){
   b_runID->GetEntry(entry);
 
   //========= initialization
-  energy = TMath::QuietNaN();
+  dE = TMath::QuietNaN();
+  TOF = TMath::QuietNaN();
+  Z = TMath::QuietNaN();
+  AoQ = TMath::QuietNaN();
   Long64_t startTimeL = 0;
   Long64_t startTimeR = 0;
-  Long64_t stopTime = 0;
+
   double cfdL = TMath::QuietNaN();
   double cfdR = TMath::QuietNaN();
-  double cfd2 = TMath::QuietNaN();
+  double cdB2 = TMath::QuietNaN();
 
+  Long64_t stopTime = 0;
   Long64_t startTime40 = 0;
   double cfd40 = TMath::QuietNaN();
 
@@ -181,40 +182,40 @@ Bool_t peachCake::Process(Long64_t entry){
     }
   }
   
-  //======= Scanning the event
+  //======= Fill detectors
   for( int i = 0; i < multi; i++){
     //----- clover
-    if( 0 <= detID[i] && detID[i] < 100 ){ 
-      if( e[i] > 0 ) {
-        int id = detID[i];
-        double eCal = 0;
-        for( int k = 0; k < int(eCorr[id].size()); k++){
-          eCal += eCorr[id][k] * TMath::Power(e[i], k);
-        }
-        hCloverID->Fill(detID[i], eCal);
-      }
-    }
-    //----- dy Beta
+    ///if( 0 <= detID[i] && detID[i] < 100 ){ 
+    ///  if( e[i] > 0 ) {
+    ///    int id = detID[i];
+    ///    double eCal = 0;
+    ///    for( int k = 0; k < int(eCorr[id].size()); k++){
+    ///      eCal += eCorr[id][k] * TMath::Power(e[i], k);
+    ///    }
+    ///    hCloverID->Fill(detID[i], eCal);
+    ///  }
+    ///}
+    //----- dy Ions (low-gain)
     if( detID[i] == 200 ) { 
       if ( multiDyBeta < 4 ){
-        dyBetaTime[multiDyBeta] = e_t[i]; 
-        dyBetaEnergy[multiDyBeta] = e[i];
-        multiDyBeta ++;
-      }
-    }
-    if( 201 <= detID[i] && detID[i] < 250){ 
-      a_H[detID[i] - 200] = e[i];
-    }
-    //----- dy Ions
-    if( detID[i] == 250 ) { 
-      if ( multiDyIons < 4 ){
-        dyIonsTime[multiDyIons] = e_t[i]; 
-        dyIonsEnergy[multiDyIons] = e[i]; 
+        dyIonsTime[multiDyBeta] = e_t[i]; 
+        dyIonsEnergy[multiDyBeta] = e[i];
         multiDyIons ++;
       }
     }
+    if( 201 <= detID[i] && detID[i] < 250){ 
+      a_L[detID[i] - 200] = e[i];
+    }
+    //----- dy Beta (high-gain)
+    if( detID[i] == 250 ) { 
+      if ( multiDyIons < 4 ){
+        dyBetaTime[multiDyIons] = e_t[i]; 
+        dyBetaEnergy[multiDyIons] = e[i]; 
+        multiDyBeta ++;
+      }
+    }
     if( 251 <= detID[i] && detID[i] < 260){
-      a_L[detID[i] - 250] = e[i];
+      a_H[detID[i] - 250] = e[i];
     }
     //----- veto_front 
     if( detID[i] == 290 ) { 
@@ -249,13 +250,13 @@ Bool_t peachCake::Process(Long64_t entry){
     //----- cross B2 logic
     if( detID[i] == 307 ) { 
       stopTime = e_t[i];
-      cfd2 = cfd[i]/16384.;
+      cdB2 = cfd[i]/16384.;
       if( (vetoFlag & 4) == 0) vetoFlag += 4;
     }
     //----- MSX40
-    ///if( detID[i] == 308 ) energy = e[i]; 
+    ///if( detID[i] == 308 ) dE = e[i]; 
     //----- MSX100
-    if( detID[i] == 309 ) energy = e[i]; 
+    if( detID[i] == 309 ) dE = e[i]; 
     //----- MSX40 logic
     if( detID[i] == 310 ) {  
       startTime40 = e_t[i];
@@ -267,38 +268,33 @@ Bool_t peachCake::Process(Long64_t entry){
   Long64_t startTime = startTimeL;
   double startCFD = cfdL;
 
-  TOF = TMath::QuietNaN();
-  if( startTime > 0 && stopTime > 0 && energy > 0 ){
+  if( startTime > 0 && stopTime > 0 && dE > 0 ){
     if( stopTime > startTime ){
-      TOF = double(stopTime - startTime) - startCFD + cfd2;
+      TOF = double(stopTime - startTime) - startCFD + cdB2;
     }else{
-      TOF = -double(startTime - stopTime) - startCFD + cfd2 ;
+      TOF = -double(startTime - stopTime) - startCFD + cdB2 ;
     }
-    TOF = TOF * 8; // to ns
+    TOF = TOF * 8; /// to ns
     
     ///============== PID correction
     if( pidCorrFileName != "" ){
       for( int k = 0; k < (int) pidCorr.size() ; k++){
         if( int(pidCorr[k][0]) != runID/100 ) continue;
-        
         TOF = pidCorr[k][1] + pidCorr[k][2] * TOF;
-        energy = pidCorr[k][3] + pidCorr[k][4] * energy + pidCorr[k][5] * energy * energy;
-        
+        dE = pidCorr[k][3] + pidCorr[k][4] * dE + pidCorr[k][5] * dE * dE;
       }
     }
     
-  
     /// TOF correction, A/Q = 3 to have a fixed TOF.
-    
-    if( energy > 160./44.*(TOF + 260) +180 ) {
-      hPID0->Fill(TOF, energy);
+    if( dE > 160./44.*(TOF + 260) +180 ) {
+      hPID0->Fill(TOF, dE);
       hTOF->Fill(TOF);
-      hdE->Fill(energy);
+      hdE->Fill(dE);
     }
 
-    double newTOF = TOF - TOFCorrection(energy) - 234;
-
-    hPID->Fill(newTOF, energy);
+    ///--------- slew correction
+    double newTOF = TOF - TOFCorrection(dE) - 234;
+    hPID->Fill(newTOF, dE);
     
     flag += 1;
 
@@ -313,18 +309,18 @@ Bool_t peachCake::Process(Long64_t entry){
     double beta = FL/c/(newTOF+tofOffset);
     double gamma = 1./TMath::Sqrt(1.-beta*beta);
 
-    double Z = sqrt((energy + 11.9473) / 23.097) * TMath::Power(beta / beta0, 1.3)  ;
-    double AoQ = c * Brho / gamma/ beta / ma;
+    Z = sqrt((dE + 11.9473) / 23.097) * TMath::Power(beta / beta0, 1.3)  ;
+    AoQ = c * Brho / gamma/ beta / ma;
 
     Z = 0.0153343 + 1.00339 * Z;
     
     hPID2->Fill(AoQ, Z);
     hZ->Fill(Z);
     
-    if(  3.5 > Z && Z >  2.5 ) hZ3->Fill( AoQ * 3);
-    if(  6.5 > Z && Z >  5.5 ) hZ6->Fill( AoQ * 6);
-    if( 10.5 > Z && Z >  9.5 ) hZ10->Fill( AoQ * 10);
-    if( 14.5 > Z && Z > 13.5 ) hZ14->Fill( AoQ * 14);
+    ///if(  3.5 > Z && Z >  2.5 ) hZ3->Fill( AoQ * 3);
+    ///if(  6.5 > Z && Z >  5.5 ) hZ6->Fill( AoQ * 6);
+    ///if( 10.5 > Z && Z >  9.5 ) hZ10->Fill( AoQ * 10);
+    ///if( 14.5 > Z && Z > 13.5 ) hZ14->Fill( AoQ * 14);
    
   }
 
@@ -341,8 +337,10 @@ Bool_t peachCake::Process(Long64_t entry){
   xBeta = (a_H[3]+a_H[2])/sumA_H;
   yBeta = (a_H[1]+a_H[2])/sumA_H;
 
+  if( 0 < xIons && xIons < 1 && 0 < yIons && yIons < 1 ) flag += 2;
+  if( 0 < xBeta && xBeta < 1 && 0 < yBeta && yBeta < 1 ) flag += 4;
 
-  if( (flag & 1) == 0 ) {
+  if( (flag & 1) == 0 ) { ///has beam
     hImplantIons->Fill(xIons, yIons);
     hImplantBeta->Fill(xBeta, yBeta);
     hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]);
@@ -353,17 +351,15 @@ Bool_t peachCake::Process(Long64_t entry){
   }
 
   if( !TMath::IsNaN(veto_f) && !TMath::IsNaN(veto_r) && crossTime == 0 ){
-    //hImplantIons->Fill(xIons, yIons);
-    //hImplantBeta->Fill(xBeta, yBeta);
-    //hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]);
-    //hImplantDyBeta->Fill(sumA_H, dyBetaEnergy[0]);
-    //
-    //hImplantX->Fill(xIons, xBeta);
-    //hImplantY->Fill(yIons, yBeta);
+    ///hImplantIons->Fill(xIons, yIons);
+    ///hImplantBeta->Fill(xBeta, yBeta);
+    ///hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]);
+    ///hImplantDyBeta->Fill(sumA_H, dyBetaEnergy[0]);
+    ///
+    ///hImplantX->Fill(xIons, xBeta);
+    ///hImplantY->Fill(yIons, yBeta);
   }
 
-  if( 0 < xIons && xIons < 1 && 0 < yIons && yIons < 1 ) flag += 2;
-  if( 0 < xBeta && xBeta < 1 && 0 < yBeta && yBeta < 1 ) flag += 4;
 
   hFlag->Fill(flag);
   hvetoFlag->Fill(vetoFlag);
@@ -375,9 +371,9 @@ Bool_t peachCake::Process(Long64_t entry){
     printf("################# %lld, multi:%d\n", entry, multi);
 
     printf("----- beam Line \n");
-    printf(" MSX100 eneergy : %f \n", energy);
+    printf(" MSX100 eneergy : %f \n", dE);
     printf(" startTime : %llu, CFD : %f \n", startTime, startCFD);
-    printf("  stopTime : %llu, CFD : %f \n", stopTime, cfd2);
+    printf("  stopTime : %llu, CFD : %f \n", stopTime, cdB2);
     printf("       TOF : %f ns \n", TOF);
 
     printf("----- Ions, low gain \n");
@@ -426,6 +422,7 @@ Bool_t peachCake::Process(Long64_t entry){
 
   //============================= Rejection
   if ( flag == 0 ) return kTRUE;
+  if ( Z < 0 ) return kTRUE;
 
   //============================= Save Tree
   if( saveNewTree ){
@@ -456,80 +453,50 @@ void peachCake::Terminate(){
     printf("|%s|\n", canvasTitle.Data());
     TCanvas * c1 = new TCanvas("c1", canvasTitle, 700 * div[0], 700 * div[1]);
     c1->Divide(div[0], div[1]);
+    int padID=0;
 
+    ///padID ++;
+    ///c1->cd(padID); //c1->cd(padID)->SetLogz();
+    ///hPID0->Draw("colz");
 
-    int padID=1;
-
-    c1->cd(padID); //c1->cd(padID)->SetLogz();
-    hPID0->Draw("colz");
-
-    //padID ++;
-    //c1->cd(padID); //c1->cd(padID)->SetLogz();
-    //hTOF->Draw();
-    //
-    //padID ++;
-    //c1->cd(padID); //c1->cd(padID)->SetLogz();
-    //hdE->Draw();
-
-    padID ++;
-    c1->cd(padID); //c1->cd(padID)->SetLogz();
-    hPID->Draw("colz");
+    ///padID ++;
+    ///c1->cd(padID); //c1->cd(padID)->SetLogz();
+    ///hPID->Draw("colz");
 
     padID ++;
     c1->cd(padID); c1->cd(padID)->SetLogz();
     hPID2->Draw("colz");
 
-    padID ++;
-    c1->cd(padID); //c1->cd(padID)->SetLogz();
-    hZ->Draw("colz");
-
-    padID ++;
-    c1->cd(padID); //c1->cd(padID)->SetLogz();
-    hZ3->Draw("colz");
-
-    padID ++;
-    c1->cd(padID); //c1->cd(padID)->SetLogz();
-    hZ6->Draw("colz");
-
-    padID ++;
-    c1->cd(padID); //c1->cd(padID)->SetLogz();
-    hZ10->Draw("colz");
-
-    padID ++;
-    c1->cd(padID); //c1->cd(padID)->SetLogz();
-    hZ14->Draw("colz");
-
     //padID ++;
     //c1->cd(padID); c1->cd(padID)->SetGrid(); c1->cd(padID)->SetLogz();
     //hCloverID->Draw("colz"); hCloverID->SetNdivisions(-13);
 
-/*
-   padID ++;
-   c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantIons->Draw("colz");
 
-   padID ++;
-   c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantBeta->Draw("colz");
-
-   //padID ++;
-   //c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantBeta_veto->Draw("colz");
-
-   padID ++;
-   c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantDyIons->Draw("colz");
-
-   padID ++;
-   c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantDyBeta->Draw("colz");
-   
-   padID ++;
-   c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantX->Draw("colz");
-   
-   padID ++;
-   c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantY->Draw("colz");
-   
-   padID ++;
-   c1->cd(padID); hFlag->Draw(""); 
-   hvetoFlag->SetLineColor(2); hvetoFlag->Draw("same");
-      
-   */
+    padID ++;
+    c1->cd(padID); c1->cd(padID)->SetLogz(); 
+    hImplantIons->Draw("colz");
+    
+    padID ++;
+    c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantBeta->Draw("colz");
+    
+    //padID ++;
+    //c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantBeta_veto->Draw("colz");
+    
+    padID ++;
+    c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantDyIons->Draw("colz");
+    
+    padID ++;
+    c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantDyBeta->Draw("colz");
+    
+    padID ++;
+    c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantX->Draw("colz");
+    
+    padID ++;
+    c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantY->Draw("colz");
+    
+    padID ++;
+    c1->cd(padID); hFlag->Draw(""); 
+    hvetoFlag->SetLineColor(2); hvetoFlag->Draw("same");
    
   }
    
diff --git a/peachCake.h b/peachCake.h
index fe320bd..7136a30 100644
--- a/peachCake.h
+++ b/peachCake.h
@@ -97,8 +97,8 @@ public :
    Int_t eventID;
    UInt_t crossEnergy;
    ULong64_t crossTime;
-   Double_t TOF;
-   Double_t energy;
+   Double_t TOF, dE;
+   Double_t Z, AoQ;
    
    Short_t flag; /// 1 = has beam, 2 = has Ions, 4 = hasBeta, default 0; 
 
@@ -161,51 +161,61 @@ void peachCake::Init(TTree *tree){
    printf(" is save tree? :%d \n", saveNewTree);
    
    if( saveNewTree ){
-     saveFileName = "test";
-     int oldSeg = -100;
-     int numFile = fChain->GetListOfFiles()->GetLast()+1;
-     for( int i = 0; i < numFile; i++){
-      TString name = fChain->GetListOfFiles()->At(i)->GetTitle();
-      int pos = name.Last('/');
-      name = name.Remove(0,pos+1); //this should be run-XXXX-XX.root
-      if( i == 0 ) {
-        saveFileName = name;
-        pos = saveFileName.Last('-');
-        saveFileName = saveFileName.Remove(pos); //run-XXXX
+      saveFileName = "zzz_";
+      int oldRun = -1;
+      bool flag1 = false;
+      int numFile = fChain->GetListOfFiles()->GetLast()+1;
+      for( int i = 0; i < numFile; i++){
+         TString name = fChain->GetListOfFiles()->At(i)->GetTitle();
+         int pos = name.Last('/');
+         name = name.Remove(0,pos+1); ///this should be run-XXXX-XX.root
+         pos = name.Last('-');
+         name = name.Remove(pos);
+         pos = name.Last('-');
+         name = name.Remove(0, pos+2);
+         int runNum = atoi(name);
+         if( runNum == oldRun + 1 ) {
+            int len = saveFileName.Sizeof();
+            if( flag1 == false) {
+               saveFileName += "-";
+            }
+            if( flag1 == true) saveFileName.Remove(len-4);
+            flag1 = true;
+         }
+         if( runNum > oldRun + 1 ) {
+            flag1 = false;
+            saveFileName += "_";
+         }
+         if( runNum > oldRun ){
+            saveFileName += Form("%03d", runNum );
+            oldRun = runNum;
+         }
       }
-      pos = name.Last('-');
-      int segNum = atoi(name.Remove(0, pos+1).Remove(2));
-      if( oldSeg == -100 ) oldSeg = segNum;
-      saveFileName = saveFileName + Form("_%02d",segNum);
-      if( segNum == oldSeg + 1 && i != numFile -1 ){
-         int len = saveFileName.Length();
-         saveFileName = saveFileName.Remove(len - 3);
-         oldSeg = segNum;
-      }  
-     }
-     saveFileName = saveFileName + ".root";
+      saveFileName = saveFileName + ".root";
      
-     printf("=========== saveFile : %s \n", saveFileName.Data());
+      printf("=========== saveFile : %s \n", saveFileName.Data());
      
-     saveFile = new TFile(saveFileName, "recreate");
-     newTree = new TTree("tree", "tree");
+      saveFile = new TFile(saveFileName, "recreate");
+      newTree = new TTree("tree", "tree");
      
-     newTree->Branch("eventID",         &eventID, "eventID/l");
-     //newTree->Branch("runID",             &runID, "runID/I");
-     newTree->Branch("TOF",                 &TOF, "TOF/D");
-     newTree->Branch("energy",           &energy, "energy/D");
-     newTree->Branch("crossTime",     &crossTime, "crossTime/l");
-     newTree->Branch("crossEnergy", &crossEnergy, "crossEnergy/i");
-     newTree->Branch("flag",               &flag, "flag/S");
-     newTree->Branch("vetoFlag",       &vetoFlag, "vetoFlag/S");
-     newTree->Branch("xIons",             &xIons, "xIons/D");
-     newTree->Branch("yIons",             &yIons, "yIons/D");
-     newTree->Branch("xBeta",             &xBeta, "xBeta/D");
-     newTree->Branch("yBeta",             &yBeta, "yBeta/D");
-     newTree->Branch("dyIonsTime",    dyIonsTime, "dyIonsTime[4]/l");
-     newTree->Branch("dyBetaTime",    dyBetaTime, "dyBetaTime[4]/l");
-     newTree->Branch("veto_f",           &veto_f, "veto_f/D");
-     newTree->Branch("veto_r",           &veto_r, "veto_r/D");     
+      newTree->Branch("eventID",         &eventID, "eventID/l");
+      newTree->Branch("runID",             &runID, "runID/I");
+      ///newTree->Branch("TOF",                 &TOF, "TOF/D");
+      ///newTree->Branch("dE",           &dE, "dE/D");
+      newTree->Branch("AoQ",              &AoQ, "AoQ/D");
+      newTree->Branch("Z",                  &Z, "Z/D");
+      newTree->Branch("crossTime",     &crossTime, "crossTime/l");
+      newTree->Branch("crossEnergy", &crossEnergy, "crossEnergy/i");
+      newTree->Branch("flag",               &flag, "flag/S");
+      newTree->Branch("vetoFlag",       &vetoFlag, "vetoFlag/S");
+      newTree->Branch("xIons",             &xIons, "xIons/D");
+      newTree->Branch("yIons",             &yIons, "yIons/D");
+      newTree->Branch("xBeta",             &xBeta, "xBeta/D");
+      newTree->Branch("yBeta",             &yBeta, "yBeta/D");
+      newTree->Branch("dyIonsTime",    dyIonsTime, "dyIonsTime[4]/l");
+      newTree->Branch("dyBetaTime",    dyBetaTime, "dyBetaTime[4]/l");
+      newTree->Branch("veto_f",           &veto_f, "veto_f/D");
+      newTree->Branch("veto_r",           &veto_r, "veto_r/D");     
    }
    
    //==== clock
diff --git a/plotPID.C b/plotPID.C
index 180258c..ecebf71 100644
--- a/plotPID.C
+++ b/plotPID.C
@@ -31,7 +31,7 @@ void plotPID() {
     
     if( nFile == 0 ) continue;
 
-    selector->SetHistRootFileName(Form("PID_%03d.root", run));
+    selector->SetHistRootFileName(Form("PID_runID/PID_%03d.root", run));
     
     chain->Process(selector, "");
     
diff --git a/script.C b/script.C
index 09d0085..39faf99 100644
--- a/script.C
+++ b/script.C
@@ -17,31 +17,31 @@ void script() {
   //chain->Add("root_data/run-0241-*.root");
   
   // new beam 
-  //chain->Add("root_data/run-024[6-8]*.root");
-  //chain->Add("root_data/run-0250*.root");
-  //chain->Add("root_data/run-025[1-4].root");
-  //chain->Add("root_data/run-025[6-9].root");
-  //chain->Add("root_data/run-026[1-4]*.root");
-  //chain->Add("root_data/run-0269*.root");
-  //chain->Add("root_data/run-027[0-1]*.root");
-  //chain->Add("root_data/run-028[3-4]*.root");
-  //chain->Add("root_data/run-028[6-8]*.root");
-  //chain->Add("root_data/run-0292*.root");
-  //chain->Add("root_data/run-029[4-6]*.root");
-  //chain->Add("root_data/run-029[8-9]*.root");
-  //chain->Add("root_data/run-031*.root");
-  //chain->Add("root_data/run-032*.root");
+  chain->Add("root_data/run-024[6-8]*.root");
+  chain->Add("root_data/run-0250*.root");
+  chain->Add("root_data/run-025[1-4].root");
+  chain->Add("root_data/run-025[6-9].root");
+  chain->Add("root_data/run-026[1-4]*.root");
+  chain->Add("root_data/run-0269*.root");
+  chain->Add("root_data/run-027[0-1]*.root");
+  chain->Add("root_data/run-028[3-4]*.root");
+  chain->Add("root_data/run-028[6-8]*.root");
+  chain->Add("root_data/run-0292*.root");
+  chain->Add("root_data/run-029[4-6]*.root");
+  chain->Add("root_data/run-029[8-9]*.root");
+  chain->Add("root_data/run-031*.root");
+  chain->Add("root_data/run-032*.root");
 
   
-  chain->Add("root_data/run-0246-00.root");
-  chain->Add("root_data/run-0246-01.root");
-  chain->Add("root_data/run-0247-00.root");
-  chain->Add("root_data/run-0247-01.root");
-  chain->Add("root_data/run-0248-00.root");
-  chain->Add("root_data/run-0250-00.root");
-  chain->Add("root_data/run-0250-01.root");
+  //chain->Add("root_data/run-0246-00.root");
+  //chain->Add("root_data/run-0246-01.root");
+  //chain->Add("root_data/run-0247-00.root");
+  //chain->Add("root_data/run-0247-01.root");
+  //chain->Add("root_data/run-0248-00.root");
+  //chain->Add("root_data/run-0250-00.root");
+  //chain->Add("root_data/run-0250-01.root");
   
-  bool isSaveNewTree = false;
+  bool isSaveNewTree = true;
   
   TString histRootFileName = "";// Form("PID_%03d.root", runNum);