Micromegas LAPP reconstruction

00001  {
00002    gROOT->Reset();
00003    gSystem->Load("/LC/Detecteurs/MicroMegas/Offline/micromegasFrameWork/trunk/lib/libMicro.so");
00004 
00005 
00006 
00007    int sigma_cut = 2; 
00008    int pedestal = 20;
00009 
00010    bool write = 1;
00011    bool print = 1;
00012 
00013 
00014 
00015 
00016    //READ INFO FROM XML
00017    std::ifstream xml;
00018    xml.open("telescope_pedestals.xml",ifstream::in);
00019    std::string line = "";
00020    std::string lineb = "";
00021    while(line.find("path=")==string::npos || lineb.find("<output")==string::npos)
00022      {
00023         lineb=line;
00024         xml>>line;
00025      }
00026    int start = line.find("\"")+1;
00027    int end =  line.rfind("\"")-6;
00028    TString file = line.substr(start,end);
00029    //TString file = "p2507_2.root";    
00030 
00031 
00032 
00033    TString output_dir = "~/micromegasFrameWork/trunk/src/analyse/TB2011/Testbeam_aug_11/Telescope/Ini_files/";
00034 
00035 
00036 
00037    TH1F * hped_mean[9];
00038    TH1F * hped_sigma[9];
00039    TH1F * hc_ped_mean[9];
00040    TH1F * hc_ped_sigma[9];
00041    TH1I * hadc[9][96];
00042    TF1  * fadc[9][96];
00043 
00044    TString name,title;
00045         
00046    for (int i=0;i<9;i++)
00047      {
00048        name = Form("hped_mean_chamber_%i",i+1);
00049        title = Form("Pedestals of chamber %i;position (ADC)",i+1);
00050        hped_mean[i] = new TH1F(name,title,200,0,400);
00051 
00052        name = Form("hc_ped_mean_chamber_%i",i+1);
00053        title = Form("Pedestals of chamber %i;channel;position (ADC)",i+1);
00054        hc_ped_mean[i] = new TH1F(name,title,96,0,96);
00055 
00056        hped_mean[i]->SetFillColor(2);
00057        hc_ped_mean[i]->SetFillColor(2);
00058 
00059        name = Form("hped_sigma_chamber_%i",i+1);
00060        title = Form("Pedestals width of chamber %i;sigma (ADC)",i+1);
00061        hped_sigma[i] = new TH1F(name,title,100,0,10);
00062 
00063        name = Form("hc_ped_sigma_chamber_%i",i+1);
00064        title = Form("Pedestals width of chamber %i;channel;sigma (ADC)",i+1);
00065        hc_ped_sigma[i] = new TH1F(name,title,96,0,96);
00066 
00067        hped_sigma[i]->SetFillColor(3);
00068        hc_ped_sigma[i]->SetFillColor(3);
00069 
00070        for (int j=0;j<96;j++)
00071          {
00072            name = Form("hadc_chamber_%i_channel_%i",i+1,j);
00073            title = Form("ADC of chamber %i channel %i",i+1,j);
00074            hadc[i][j] = new TH1I(name,title,1024,0,1024);
00075            hadc[i][j]->SetFillColor(4);
00076 
00077            name = Form("fadc_chamber_%i_channel_%i",i+1,j);
00078            fadc[i][j] = new TF1(name,"gaus",0,1024);
00079            fadc[i][j]->SetLineColor(1);
00080            fadc[i][j]->SetLineWidth(2);
00081          }
00082      }
00083 
00084 
00085 
00086 
00087 
00088 
00089    /*READ ROOT FILE*/
00090 
00091 
00092    cout<<""<<endl;
00093    cout<<"READ DATA FROM FILE "<<file<<endl;
00094    cout<<""<<endl;
00095 
00096 
00097    TFile *f = new TFile(file);
00098 
00099    TIter nextkey(f.GetListOfKeys());  
00100    TKey *key;
00101         
00102    int nchannel = 0;
00103    int adc = 0;
00104    int hardid = 0;
00105    int chbid = 0;       
00106 
00107    while (key = (TKey*)nextkey()) 
00108      {
00109     
00110        TTree *tree = (TTree*)key->ReadObj();                            
00111        MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
00112 
00113        int nEvent = tree.GetEntries();   
00114        cout <<"Number of events = "<<nEvent<<endl;
00115        cout<<endl;
00116 
00117        MTEvent *evt = new MTEvent();
00118        TBranch *branch= tree->GetBranch("MTEvent");
00119        branch->SetAddress(&evt);
00120 
00121        for (int i=1;i<nEvent+1;i++)
00122          {
00123 
00124            tree.GetEntry(i);
00125            nchannel = evt->GetNchannel();
00126            
00127            for (int j=0;j<nchannel;j++)
00128              { 
00129 
00130                MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00131                chbid = channel->GetChamberId();
00132                hardid = channel->GetHardId();
00133                adc = channel->GetAnalogValue();
00134                hadc[chbid-1][hardid]->Fill(adc);
00135 
00136              }
00137          }
00138      }
00139 
00140         
00141 
00142 
00143    /*FIT PEDESTAL DISTRIBUTIONS*/
00144 
00145    cout<<""<<endl;
00146    cout<<"FIT PEDESTAL DISRIBUTIONS"<<endl;
00147    cout<<""<<endl;
00148 
00149    for (int i=0;i<9;i++)
00150      {
00151        for (int j=0;j<96;j++)
00152          {
00153            if (hadc[i][j]->GetEntries()!=0)
00154              {
00155 
00156                fadc[i][j]->SetParameters(hadc[i][j]->GetMaximum(),hadc[i][j]->GetMean(),hadc[i][j]->GetRMS());
00157                hadc[i][j]->Fit(fadc[i][j],"q");
00158 
00159                hped_mean[i]->Fill(fadc[i][j]->GetParameter(1));
00160                hc_ped_mean[i]->SetBinContent(j+1,fadc[i][j]->GetParameter(1));
00161                hc_ped_mean[i]->SetBinError(j+1,fadc[i][j]->GetParError(1));
00162 
00163                hped_sigma[i]->Fill(fadc[i][j]->GetParameter(2));
00164                hc_ped_sigma[i]->SetBinContent(j+1,fadc[i][j]->GetParameter(2));
00165                hc_ped_sigma[i]->SetBinError(j+1,fadc[i][j]->GetParError(2));
00166 
00167              }
00168          }
00169      }
00170 
00171 
00172 
00173 
00174    /*WRITE TO FILE*/
00175 
00176    if (write){
00177 
00178      ofstream out1("Ini_files/PedestalADC1.ini");
00179      ofstream out2("Ini_files/PedestalADC2.ini");
00180      ofstream out3("Ini_files/PedestalADC3.ini");
00181      ofstream out4("Ini_files/PedestalADC4.ini");
00182      ofstream out5("Ini_files/PedestalADC5.ini");
00183      
00184      
00185      ofstream outcut1("Ini_files/PedestalADC1_cut.ini");
00186      ofstream outcut2("Ini_files/PedestalADC2_cut.ini");
00187      ofstream outcut3("Ini_files/PedestalADC3_cut.ini");
00188      ofstream outcut4("Ini_files/PedestalADC4_cut.ini");
00189      ofstream outcut5("Ini_files/PedestalADC5_cut.ini");
00190 
00191         
00192      /*First ADC module*/
00193 
00194      i = 0;
00195      out1<<"BLOCK0"<<endl;
00196      outcut1<<"BLOCK0"<<endl;
00197      for (int j=0;j<96;j++)
00198        {
00199          out1<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<0<<endl;
00200          outcut1<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<floor(0.5+fadc[i][j]->GetParameter(1) - pedestal + fadc[i][j]->GetParameter(2) * sigma_cut)<<endl;
00201        }
00202      out1<<"-1"<<endl;
00203      outcut1<<"-1"<<endl;
00204      
00205      i = 1;
00206      out1<<"BLOCK1"<<endl;
00207      outcut1<<"BLOCK1"<<endl;
00208      for (int j=0;j<96;j++)
00209        {
00210          out1<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<0<<endl;
00211          outcut1<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<floor(0.5+fadc[i][j]->GetParameter(1) - pedestal + fadc[i][j]->GetParameter(2) * sigma_cut)<<endl;
00212        }
00213      out1<<"-1"<<endl;
00214      outcut1<<"-1"<<endl;
00215      
00216      
00217      
00218      /*Second ADC module*/
00219      i = 2;
00220      out2<<"BLOCK0"<<endl;
00221      outcut2<<"BLOCK0"<<endl;
00222      for (int j=0;j<96;j++)
00223        {
00224          out2<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<0<<endl;
00225          outcut2<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<floor(0.5+fadc[i][j]->GetParameter(1) - pedestal + fadc[i][j]->GetParameter(2) * sigma_cut)<<endl;
00226        }
00227      out2<<"-1"<<endl;
00228      outcut2<<"-1"<<endl;
00229      
00230      i = 3;
00231      out2<<"BLOCK1"<<endl;
00232      outcut2<<"BLOCK1"<<endl;
00233      for (int j=0;j<96;j++)
00234        {
00235          out2<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<0<<endl;
00236          outcut2<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<floor(0.5+fadc[i][j]->GetParameter(1) - pedestal + fadc[i][j]->GetParameter(2) * sigma_cut)<<endl;
00237        }
00238      out2<<"-1"<<endl;
00239      outcut2<<"-1"<<endl;
00240 
00241 
00242 
00243      /*Third ADC module*/
00244      i = 4;
00245      out3<<"BLOCK0"<<endl;
00246      outcut3<<"BLOCK0"<<endl;
00247      for (int j=0;j<96;j++)
00248        {
00249          out3<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<0<<endl;
00250          outcut3<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<floor(0.5+fadc[i][j]->GetParameter(1) - pedestal + fadc[i][j]->GetParameter(2) * sigma_cut)<<endl;
00251        }
00252      out3<<"-1"<<endl;
00253      outcut3<<"-1"<<endl;
00254      
00255      i = 5;
00256      out3<<"BLOCK1"<<endl;
00257      outcut3<<"BLOCK1"<<endl;
00258      for (int j=0;j<96;j++)
00259        {
00260          out3<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<0<<endl;
00261          outcut3<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<floor(0.5+fadc[i][j]->GetParameter(1) - pedestal + fadc[i][j]->GetParameter(2) * sigma_cut)<<endl;
00262        }
00263      out3<<"-1"<<endl;
00264      outcut3<<"-1"<<endl;
00265      
00266      
00267 
00268      /*Fourth ADC module*/
00269      i = 6;
00270      out4<<"BLOCK0"<<endl;
00271      outcut4<<"BLOCK0"<<endl;
00272      for (int j=0;j<96;j++)
00273        {
00274          out4<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<0<<endl;
00275          outcut4<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<floor(0.5+fadc[i][j]->GetParameter(1) - pedestal + fadc[i][j]->GetParameter(2) * sigma_cut)<<endl;
00276        }
00277      out4<<"-1"<<endl;
00278      outcut4<<"-1"<<endl;
00279      
00280      i = 7;
00281      out4<<"BLOCK1"<<endl;
00282      outcut4<<"BLOCK1"<<endl;
00283      for (int j=0;j<96;j++)
00284        {
00285          out4<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<0<<endl;
00286          outcut4<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<floor(0.5+fadc[i][j]->GetParameter(1) - pedestal + fadc[i][j]->GetParameter(2) * sigma_cut)<<endl;
00287        }
00288      out4<<"-1"<<endl;
00289      outcut4<<"-1"<<endl;
00290 
00291 
00292 
00293      /*Fith ADC module*/
00294      i = 8;
00295      out5<<"BLOCK0"<<endl;
00296      outcut5<<"BLOCK0"<<endl;
00297      for (int j=0;j<96;j++)
00298        {
00299          out5<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<0<<endl;
00300          outcut5<<j<<","<<floor(0.5+fadc[i][j]->GetParameter(1)) - pedestal<<","<<floor(0.5+fadc[i][j]->GetParameter(1) - pedestal + fadc[i][j]->GetParameter(2) * sigma_cut)<<endl;
00301        }
00302      out5<<"-1"<<endl;
00303      outcut5<<"-1"<<endl;
00304 
00305    }
00306 
00307 
00308 
00309 
00310    if (print)
00311      {
00312 
00313        cout<<endl;
00314        cout<<"PRINT PLOTS"<<endl;
00315        cout<<endl;
00316 
00317 
00318        /*Strip chambers 1-2*/
00319 
00320        TCanvas * c1 = new TCanvas("c1","Pedestal alignment telescope",10,10,1000,1000);
00321        c1->Divide(2,2);
00322 
00323        for (int i=0;i<2;i++)
00324          {
00325            c1->cd(i+1);
00326            hped_mean[i]->Draw();
00327            c1->cd(i+1+2);
00328            hc_ped_mean[i]->Draw();
00329          }
00330 
00331        c1->Print("Plots/hped_mean_strip_chb_12.png");
00332 
00333        for (int i=0;i<2;i++)
00334          {
00335            c1->cd(i+1);
00336            hped_sigma[i]->Draw();
00337            c1->cd(i+1+2);
00338            hc_ped_sigma[i]->Draw();
00339          }
00340 
00341        c1->Print("Plots/hped_sigma_strip_chb_12.png");
00342 
00343 
00344        /*Strip chambers 3-4*/
00345 
00346        for (int i=0;i<2;i++)
00347          {
00348            c1->cd(i+1);
00349            hped_mean[i+2]->Draw();
00350            c1->cd(i+1+2);
00351            hc_ped_mean[i+2]->Draw();
00352          }
00353 
00354        c1->Print("Plots/hped_mean_strip_chb_34.png");
00355 
00356        for (int i=0;i<2;i++)
00357          {
00358            c1->cd(i+1);
00359            hped_sigma[i+2]->Draw();
00360            c1->cd(i+1+2);
00361            hc_ped_sigma[i+2]->Draw();
00362          }
00363 
00364        c1->Print("Plots/hped_sigma_strip_chb_34.png");
00365 
00366 
00367        /*Strip chambers 5-6*/
00368 
00369        TCanvas * c1 = new TCanvas("c1","Pedestal alignment telescope",10,10,1000,1000);
00370        c1->Divide(2,2);
00371 
00372        for (int i=0;i<2;i++)
00373          {
00374            c1->cd(i+1);
00375            hped_mean[i+4]->Draw();
00376            c1->cd(i+1+2);
00377            hc_ped_mean[i+4]->Draw();
00378          }
00379 
00380        c1->Print("Plots/hped_mean_strip_chb_56.png");
00381 
00382        for (int i=0;i<2;i++)
00383          {
00384            c1->cd(i+1);
00385            hped_sigma[i+4]->Draw();
00386            c1->cd(i+1+2);
00387            hc_ped_sigma[i+4]->Draw();
00388          }
00389 
00390        c1->Print("Plots/hped_sigma_strip_chb_56.png");
00391 
00392 
00393 
00394 
00395        /*Pad chambers*/
00396 
00397        TCanvas * c0 = new TCanvas("c0","Pedestal alignment telescope",10,10,1500,1000);
00398        c0->Divide(3,2);
00399 
00400        for (int i=0;i<3;i++)
00401          {
00402            c0->cd(i+1);
00403            hped_mean[i+6]->Draw();
00404            c0->cd(i+1+3);
00405            hc_ped_mean[i+6]->Draw();
00406          }
00407 
00408        c0->Print("Plots/hped_mean_pad_chb_789.png");
00409 
00410        for (int i=0;i<3;i++)
00411          {
00412            c0->cd(i+1);
00413            hped_sigma[i+6]->Draw();
00414            c0->cd(i+1+3);
00415            hc_ped_sigma[i+6]->Draw();
00416          }
00417 
00418        c0->Print("Plots/hped_sigma_pad_chb_789.png");
00419 
00420 
00421 
00422 
00423 
00424 
00425 
00426      }
00427 
00428    
00429 
00430  }
00431 
00432 
00433
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/TB2011/Testbeam_aug_11/Telescope/Pedestals/telescope_pedestals.C
