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