/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/Analysis_dirac.C

Go to the documentation of this file.

{gROOT->Reset();
  
  gSystem->Load("libMicro.so");
  gStyle->SetPalette(1);

  bool display = 0;
  bool plot = 1;

  TString file = "dirac";

  TString filename =  file +".root";

  int xpos,ypos,zpos,val;
  Long64_t t,t1,t2,t3,dt;
  int nchannel = 0;
  bool time_cut = 0;
  int hardid = 0;
  int softid = 0;
  int nhit = 0;

  TH1I * hx_1 = new TH1I("hx_1","X distribution of hits in chamber 1",8,0,8);
  TH1I * hx_2 = new TH1I("hx_2","X distribution of hits in chamber 2",8,0,8);
  TH1I * hx_3 = new TH1I("hx_3","X distribution of hits in chamber 3",8,0,8);
  TH1I * hx_4 = new TH1I("hx_4","X distribution of hits in chamber 4",8,0,8);

  TH1I * hy_1 = new TH1I("hy_1","Y distribution of hits in chamber 1",8,0,8);
  TH1I * hy_2 = new TH1I("hy_2","Y distribution of hits in chamber 2",8,0,8);
  TH1I * hy_3 = new TH1I("hy_3","Y distribution of hits in chamber 3",8,0,8);
  TH1I * hy_4 = new TH1I("hy_4","Y distribution of hits in chamber 4",8,0,8);

  TH2I * hxy_1 = new TH2I("hxy_1","XY distribution of hits in chamber 1",8,0,8,8,0,8);
  TH2I * hxy_2 = new TH2I("hxy_2","XY distribution of hits in chamber 2",8,0,8,8,0,8);
  TH2I * hxy_3 = new TH2I("hxy_3","XY distribution of hits in chamber 3",8,0,8,8,0,8);
  TH2I * hxy_4 = new TH2I("hxy_4","XY distribution of hits in chamber 4",8,0,8,8,0,8);

  TH1I * hz = new TH1I("hz","Z distribution of hits",25,-5,20);

  TH1I * hc_1 = new TH1I("hc_1","Channel number of hits in chamber 1",64,0,64);
  TH1I * hc_2 = new TH1I("hc_2","Channel number of hits in chamber 2",64,0,64);
  TH1I * hc_3 = new TH1I("hc_3","Channel number of hits in chamber 3",64,0,64);
  TH1I * hc_4 = new TH1I("hc_4","Channel number of hits in chamber 4",64,0,64);

  int bcid_hit_max = 10000;
  TH1I * hbcid_hit_1 = new TH1I("hbcid_hit_1","BCID_hit distribution in chamber 1",1000,0,bcid_hit_max);
  TH1I * hbcid_hit_2 = new TH1I("hbcid_hit_2","BCID_hit distribution in chamber 2",1000,0,bcid_hit_max);
  TH1I * hbcid_hit_3 = new TH1I("hbcid_hit_3","BCID_hit distribution in chamber 3",1000,0,bcid_hit_max);
  TH1I * hbcid_hit_4 = new TH1I("hbcid_hit_4","BCID_hit distribution in chamber 4",1000,0,bcid_hit_max);

  int bcid_dif_max = 2000;
  TH1I * hbcid_dif_1 = new TH1I("hbcid_dif_1","BCID_dif distribution in chamber 1",1000,0,bcid_dif_max);
  TH1I * hbcid_dif_2 = new TH1I("hbcid_dif_2","BCID_dif distribution in chamber 2",1000,0,bcid_dif_max);
  TH1I * hbcid_dif_3 = new TH1I("hbcid_dif_3","BCID_dif distribution in chamber 3",1000,0,bcid_dif_max);
  TH1I * hbcid_dif_4 = new TH1I("hbcid_dif_4","BCID_dif distribution in chamber 4",1000,0,bcid_dif_max);

  int bcid_abs_max = 1e9;
  TH1I * hbcid_abs_1 = new TH1I("hbcid_abs_1","BCID_abs distribution in chamber 1",1000,0,bcid_abs_max);
  TH1I * hbcid_abs_2 = new TH1I("hbcid_abs_2","BCID_abs distribution in chamber 2",1000,0,bcid_abs_max);
  TH1I * hbcid_abs_3 = new TH1I("hbcid_abs_3","BCID_abs distribution in chamber 3",1000,0,bcid_abs_max);
  TH1I * hbcid_abs_4 = new TH1I("hbcid_abs_4","BCID_abs distribution in chamber 4",1000,0,bcid_abs_max);

  TH1I * hnhit_1 = new TH1I("hnhit_1","Number of hits per event (per DIF readout) in chamber 1",100,0,100);
  TH1I * hnhit_2 = new TH1I("hnhit_2","Number of hits per event (per DIF readout) in chamber 2",100,0,100);
  TH1I * hnhit_3 = new TH1I("hnhit_3","Number of hits per event (per DIF readout) in chamber 3",100,0,100);
  TH1I * hnhit_4 = new TH1I("hnhit_4","Number of hits per event (per DIF readout) in chamber 4",100,0,100);



  int n1=0;
  int n2=0;
  int n3=0;
  int n4=0;
  int one_hit_each_chamber = 0;



  TFile *f = new TFile(filename);
  cout<<"reading File : "<<filename<<endl;


  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;
  float nevent = nEvent;


  MTEvent *evt = new MTEvent();

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





  for (int i=0;i<nevent;i++){  
    //for (int i=nevent-10;i<nevent;i++){  

    if (i%10 == 0){cout<<i/nevent*100<<" %"<<"\r"<<flush;}

    tree.GetEntry(i);
    nchannel = evt->GetNchannel();
    nhit+=nchannel;

    //cout<<"event number = "<<i<<"  Number of hits = "<<nchannel<<endl;


    n1=0;n2=0;n3=0;n4=0;


    for (int j=0;j<nchannel;j++){

    MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
    xpos = channel->GetX();
    ypos = channel->GetY();
    zpos = channel->GetZ();
    val = channel->GetDigitalValue();
    hardid = channel->GetHardId();
    softid = channel->GetSoftId();

    t1 = channel->GetBcId_Abs();
    t2 = channel->GetBcId_Dif();
    t3 = channel->GetBcId_Hit();
    dt = (t2-t3);
    t = (t1 - (t2-t3));


    hz->Fill(zpos);

    if (zpos == 0) {
      hx_1->Fill(xpos);hy_1->Fill(ypos);hxy_1->Fill(xpos,ypos);hc_1->Fill(hardid);hbcid_hit_1->Fill(t3);hbcid_dif_1->Fill(t2);hbcid_abs_1->Fill(t1);hnhit_1->Fill(nchannel);if (hardid!=11 && hardid!=57){n1++;}}

    if (zpos == 5) {
      hx_2->Fill(xpos);hy_2->Fill(ypos);hxy_2->Fill(xpos,ypos);hc_2->Fill(hardid);hbcid_hit_2->Fill(t3);hbcid_dif_2->Fill(t2);hbcid_abs_2->Fill(t1);hnhit_2->Fill(nchannel);if (hardid!=11 && hardid!=57){n2++;}}

    if (zpos == 10){
      hx_3->Fill(xpos);hy_3->Fill(ypos);hxy_3->Fill(xpos,ypos);hc_3->Fill(hardid);hbcid_hit_3->Fill(t3);hbcid_dif_3->Fill(t2);hbcid_abs_3->Fill(t1);hnhit_3->Fill(nchannel);if (hardid!=11 && hardid!=57){n3++;}}

    if (zpos == 15){
      hx_4->Fill(xpos);hy_4->Fill(ypos);hxy_4->Fill(xpos,ypos);hc_4->Fill(hardid);hbcid_hit_4->Fill(t3);hbcid_dif_4->Fill(t2);hbcid_abs_4->Fill(t1);hnhit_4->Fill(nchannel);if (hardid!=11 && hardid!=57){n4++;}}

    time_cut = (dt>=2 && dt<=8);

    if (display){cout<<" x = "<<xpos<<" y = "<<ypos<<" z = "<<zpos<<" val = "<<val<<" hardid = "<<hardid<<" softid = "<<softid<<" bcid_abs = "<<t1<<" bcid_dif = "<<t2<<" bcid_hit = "<<t3<<endl;}}}}

   


  cout<<"Number of hits = "<<nhit<<endl;




  if (plot){

    TCanvas * c1 = new TCanvas("c1","",10,10,1600,400);
    c1->Divide(4,1);
    c1->cd(1);
    hxy_1->Draw("zcol");
    c1->cd(2);
    hxy_2->Draw("zcol");
    c1->cd(3);
    hxy_3->Draw("zcol");
    c1->cd(4);
    hxy_4->Draw("zcol");

    TCanvas * c2 = new TCanvas("c2","",10,10,1600,400);
    c2->Divide(4,1);
    c2->cd(1);
    hc_1->Draw();
    c2->cd(2);
    hc_2->Draw();
    c2->cd(3);
    hc_3->Draw();
    c2->cd(4);
    hc_4->Draw();

    TCanvas * c3 = new TCanvas("c3","",10,10,1600,400);
    c3->Divide(4,1);
    c3->cd(1);
    hbcid_hit_1->Draw();
    c3->cd(2);
    hbcid_hit_2->Draw();
    c3->cd(3);
    hbcid_hit_3->Draw();
    c3->cd(4);
    hbcid_hit_4->Draw();

    TCanvas * c4 = new TCanvas("c4","",10,10,600,600);
    hz->Draw();

    TCanvas * c5 = new TCanvas("c5","",10,10,1600,400);
    c5->Divide(4,1);
    c5->cd(1);
    hnhit_1->Draw();
    c5->cd(2);
    hnhit_2->Draw();
    c5->cd(3);
    hnhit_3->Draw();
    c5->cd(4);
    hnhit_4->Draw();

    TCanvas * c6 = new TCanvas("c6","",10,10,1600,400);
    c6->Divide(4,1);
    c6->cd(1);
    hbcid_dif_1->Draw();
    c6->cd(2);
    hbcid_dif_2->Draw();
    c6->cd(3);
    hbcid_dif_3->Draw();
    c6->cd(4);
    hbcid_dif_4->Draw();

    TCanvas * c7 = new TCanvas("c7","",10,10,1600,400);
    c7->Divide(4,1);
    c7->cd(1);
    hbcid_abs_1->Draw();
    c7->cd(2);
    hbcid_abs_2->Draw();
    c7->cd(3);
    hbcid_abs_3->Draw();
    c7->cd(4);
    hbcid_abs_4->Draw();


    TString plotname1 = "Plots/"+ file +"_xy_hit.png";
    TString plotname2 = "Plots/"+ file +"_channel_number_hit.png";
    TString plotname3 = "Plots/"+ file +"_bcid_hit.png";
    TString plotname4 = "Plots/"+ file +"_z_hit.png";
    TString plotname5 = "Plots/"+ file +"_number_of_hit.png";
    TString plotname6 = "Plots/"+ file +"_bcid_dif.png";
    TString plotname7 = "Plots/"+ file +"_bcid_abs.png";

    c1->Print(plotname1.Data());
    c2->Print(plotname2.Data());
    c3->Print(plotname3.Data());
    c4->Print(plotname4.Data());
    c5->Print(plotname5.Data());
    c6->Print(plotname6.Data());
    c7->Print(plotname7.Data());



  }




}

---