musrsim/trimsp/TrimSPGUI/plotRge.C

// File:        plotRge.C
// Author:      Zaher Salman, originally by Andreas Suter
// Date:        03/06/2009
// Purpose:     ROOT macro to read and plot Range data from trimsp calculation
//              Assume following file Format:
//
//  $Id$
//
//  DEPTH      PARTICLES
//    15.         260
//    45.         595
//    75.         820
//   105.        1179
//
//  plotRge(names), e.g.
//  plotRge("InSne_E2500, InSne_E4000");
//  will read files

void plotRge(TString &names)
{
  if (names.CompareTo("help", TString::kIgnoreCase) == 0) {
    cout << endl << "usage: plotRge(<names>) or plotRge(\"help\")";
    cout << endl << "  <names> is a list of rge-files";
    cout << endl << endl;
    return;
  }

  TObjArray *tokens;
  TObjString *ostr;
  TString str;
  char cstr[128];

  Int_t i, j;
  Int_t nStep = 0;
  Double_t *depth;
  Double_t *nStop;
  Double_t norm;
  Double_t intNorm;
  Double_t normStop[1000];

  TGraph *rge[1000];

  tokens = names.Tokenize(" \t,;");

  Double_t min = 0, max = 0, xmax = 0, xmin = 0;
  Int_t imax;

  // Loop over file names
  for (i=0; i<tokens->GetEntries(); i++) {
    ostr = dynamic_cast<TObjString*>(tokens->At(i));
    str = ostr->GetString();
    cout << endl << i << ": read " << str.Data();

    TGraph *fileGraph = new TGraph(str);
    if ( fileGraph->IsZombie() ){
      cout << endl;
      cout << "File " << str << " does not exist!!!" << endl << endl;
      return;
    }

    depth = fileGraph->GetX();
    nStop = fileGraph->GetY();
    nStep = fileGraph->GetN(); // Number of points in range

    // calculate step size in nm (convert from AA)
    norm = (depth[1] - depth[0])/10;

    // intNorm is the total stopping particles
    intNorm = 0.0;
    for (j=0; j<nStep; j++) {
      intNorm += nStop[j];
    }


    for (j=0; j<nStep; j++) {
      // Calculate normalized stopping profile
      if (norm > 0.) {
	// devide original profile by (step size * total particles) in %
        normStop[j] =  100*nStop[j]/(norm*intNorm);
      } else {
        normStop[j] = -1.;
      }
      // rescale depth to nm
      depth[j] = depth[j]/10.;  

      // Find max and min of stopping profiles
      if (normStop[j] > max)
        max = normStop[j];
      if (normStop[j] < min)
        min = normStop[j];

      // Find max and min of stopping depth
      if (depth[j] > xmax) {
        xmax = depth[j];
	imax = i;
      }
      if (depth[j] < xmin) {
        xmin = depth[j];
      }
    }

    rge[i] = new TGraph(nStep, depth, normStop);

    ostr = dynamic_cast<TObjString*>(tokens->At(i));
    str  = ostr->GetString() + ", ";
    sprintf(cstr, "%.1lf", ostr->GetString().Atof() /1000.0);
    str += cstr;
    str += " (keV)";
    rge[i]->SetTitle(str.Data());
    rge[i]->SetFillColor(TColor::kWhite);

    delete fileGraph;
  }
  cout << endl;
  // plot again deepest to get the right x scale !!
  // rge[imax] = new TGraph(nStep, depth, normStop);

  Int_t color[1000];
  for (i=0; i<1000; i++)
    color[i] = -1;
  color[0] = TColor::kRed;
  color[1] = TColor::kGreen;
  color[2] = TColor::kBlue;
  color[3] = TColor::kMagenta;
  color[4] = TColor::kOrange;
  color[5] = TColor::kViolet;
  color[6] = TColor::kAzure+7;
  color[7] = TColor::kOrange+4;
  color[8] = TColor::kBlue-7;

  Int_t last = tokens->GetEntries()-1;
  // Int_t last = imax;
  rge[last]->SetMarkerStyle(20);
  rge[last]->GetXaxis()->SetTitle("Depth [nm]");
  rge[last]->GetYaxis()->SetTitle("Normalized Stopping Profile [%/nm]");

  TCanvas *c1 = new TCanvas("c1", "Range");
  c1->Show();
  
  rge[last]->Draw("apc");
  rge[last]->GetYaxis()->SetRangeUser(min, 1.05*max);
  rge[last]->GetXaxis()->SetRangeUser(xmin, xmax);
    //  cout << "xmin " << xmin << " xmax" << xmax << endl;


  Int_t col;
  for (i=0; i<last; i++) {
    if (color[i] >= 0) {
      rge[i]->SetMarkerColor(color[i]);
      rge[i]->SetLineColor(color[i]);
    } else {
      TRandom *rand = new TRandom(i);
      col = TColor::GetColor((Int_t)rand->Integer(255),(Int_t)rand->Integer(255),(Int_t)rand->Integer(255));
      rge[i]->SetMarkerColor(col);
      rge[i]->SetLineColor(col);
      delete rand;
    }
    rge[i]->SetMarkerStyle(20);
    rge[i]->Draw("plsame");
  }

  TLegend *legend = new TLegend(0.7, 0.7, 0.95, 0.95);
  legend->SetFillColor(TColor::GetColor(255,255,255)); // white
  for (i=0; i<tokens->GetEntries(); i++) {
    legend->AddEntry(rge[i]);
  }
  legend->Draw();

  delete tokens;

  // I am not sure what this does, but it waits until canvas is closed
  c1->WaitPrimitive(" ");
  cout << endl << "Canvas Closed" << endl ;

  // Then quit root cleanly
  gApplication->Terminate();

}