/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/TB2010/Testbeam_june_10/Gassiplex/Beam_profile/Beam_profile.C

Go to the documentation of this file.

{gROOT->Reset();
  
  gSystem->Load("/LC/Detecteurs/MicroMegas/Offline/micromegasFrameWork/trunk/libMicro.so");
  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);


  bool plot = 1;
  bool fit = 1;
  


  int cut = 30;




  //read file info from xml
  std::ifstream xml;
  xml.open("gassiplex_beam_profile.xml",ifstream::in);
  std::string line="";
  std::string lineb="";
  while(line.find("path=")==string::npos || lineb.find("<output")==string::npos)
   {
        lineb=line;
        xml>>line;
   }
  int start = line.find("\"")+1;
  int end =  line.rfind("\"")-6;
  TString address = line.substr(start,end);




  int nchannel = 0;
  float xpos = 0;
  float ypos = 0;
  float adc = 0;
  int content = 0;

  TH2F ** hxy = new TH2F[4];
  for(int chb = 0;chb<4;chb++){hxy[chb] = new TH2F(Form("hxy_%i",chb),Form("Chamber %i;y (cm);x (cm)",chb),32,0,32,12,0,12);}
  




  cout << "open rootfile "<<address<<endl;
  TFile *f = new TFile(address);

  TIter nextkey(f.GetListOfKeys());
  TKey *key;


  while (key = (TKey*)nextkey()) {
    
    TTree *tree = (TTree*)key->ReadObj();                
    MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
    //run->Info();

    int nEvent = tree.GetEntries();
    cout <<"nEvent="<<nEvent<<endl;


  MTEvent *evt = new MTEvent();

  TBranch *branch= tree->GetBranch("MTEvent");
  branch->SetAddress(&evt);


  int chamber;
  int progress=0;

  
  for (int i=0;i<nEvent;i++){  
    
    if(progress!=(100*i)/nEvent){progress=(100*i/nEvent);cout<<"Progress "<<progress+1<<" %\t\r"<<flush;}
    tree.GetEntry(i);
    nchannel = evt->GetNchannel();
    
    for(int chan=0;chan<nchannel;chan++){
      
      MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(chan);
      
      chamber = channel->GetChamberId()-1;
      adc = channel->GetAnalogValue();
      xpos = channel->GetX();
      ypos = channel->GetY();
      
      if(adc>cut){hxy[chamber]->Fill(ypos,xpos);}}}}
 





  if (plot){

    TCanvas * c0 = new TCanvas("c0","",10,10,700,700);
    c0->Divide(2,2);

    for(int chb = 0;chb<4;chb++){

        c0->cd(chb+1);
        hxy[chb]->Draw("zcol");}}





if (fit){


 hx1 = hxy[0]->ProjectionX();
 hy1 = hxy[0]->ProjectionY();

 hx2 = hxy[1]->ProjectionX();
 hy2 = hxy[1]->ProjectionY();

 hx3 = hxy[2]->ProjectionX();
 hy3 = hxy[2]->ProjectionY();

 hx4 = hxy[3]->ProjectionX();
 hy4 = hxy[3]->ProjectionY();}}

---