Micromegas LAPP reconstruction

00001 {gROOT->Reset();
00002   
00003   gSystem->Load("/LC/framework_tb2010/libMicro.so");
00004   
00005 
00006 
00007 //  bool cut = 1;
00008   TString sigma_cut = "2";
00009 
00010   int pedestal = 20;
00011 
00012 
00013   TString output_path = "/LC/framework_tb2010/pedestal_gassiplex/";
00014 
00015   //TString address = "/lapp_data/LC/Detecteurs/MicroMegas/data/TB2010/SPS_H4_june_2010/Root_files/.root";
00016 
00017   bool display = 0;
00018 
00019 
00020   std::ifstream xml;
00021   xml.open("gassiplex_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   TString address = line.substr(6,line.find("\"/>")-6);
00030 
00031   cout<<"reading from "<<address<<endl;
00032 
00033 
00034 
00035 
00036   TH1I * hchamber_0[96];
00037   TH1I * hchamber_1[96];
00038   TH1I * hchamber_2[96];
00039   TH1I * hchamber_3[96*4];
00040 
00041   TString name,title;
00042 
00043   for (int i=0;i<96;i++){
00044 
00045     name = Form("hchamber_0_channel_%i",i);
00046     title = Form("ADC of chamber_0_channel_%i",i);
00047     hchamber_0[i] = new TH1I(name,title,1024,0,1024);
00048 
00049     name = Form("hchamber_1_channel_%i",i);
00050     title = Form("ADC of chamber_1_channel_%i",i);
00051     hchamber_1[i] = new TH1I(name,title,1024,0,1024);
00052 
00053     name = Form("hchamber_2_channel_%i",i);
00054     title = Form("ADC of chamber_2_channel_%i",i);
00055     hchamber_2[i] = new TH1I(name,title,1024,0,1024);}
00056 
00057   for (int i=0;i<96*4;i++){
00058 
00059     name = Form("hchamber_3_channel_%i",i);
00060     title = Form("ADC of chamber_3_channel_%i",i);
00061     hchamber_3[i] = new TH1I(name,title,1024,0,1024);}
00062 
00063 
00064 cout << "open rootfile"<<endl;
00065   TFile *f = new TFile(address);
00066 
00067   TIter nextkey(f.GetListOfKeys());
00068   TKey *key;
00069 
00070   while (key = (TKey*)nextkey()) {
00071     
00072     TTree *tree = (TTree*)key->ReadObj();                
00073     MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
00074     //run->Info();
00075 
00076   int nEvent = tree.GetEntries();
00077   cout <<"nEvent="<<nEvent<<endl;
00078 
00079 
00080   MTEvent *evt = new MTEvent();
00081 
00082   TBranch *branch= tree->GetBranch("MTEvent");
00083   branch->SetAddress(&evt);
00084 
00085   int nchannel = 0;
00086   float zpos = 0;
00087   float adc = 0;
00088   int hardid = 0;
00089 
00090 
00091   cout<<endl;
00092   cout<<"Reading"<<endl;
00093 
00094   int progress = 0;
00095   for (int i=1;i<nEvent;i++){  
00096   if(progress!=100*i/nEvent){
00097         progress=100*i/nEvent;
00098         cout<<"progress "<<progress+1<<" % \t\r"<<flush;
00099         }
00100   
00101     tree.GetEntry(i);
00102     nchannel = evt->GetNchannel();
00103 
00104     cout<<"event number = "<<i<<"  Number of hits = "<<nchannel<<endl;
00105 
00106     for (int j=0;j<96;j++){
00107       MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00108       zpos = channel->GetZ();
00109       hardid = channel->GetHardId();
00110       adc = channel->GetAnalogValue();
00111       hchamber_0[hardid]->Fill(adc);}
00112 
00113     for (int j=96;j<96*2;j++){
00114       MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00115       zpos = channel->GetZ();
00116       hardid = channel->GetHardId();
00117       adc = channel->GetAnalogValue();
00118       hchamber_1[hardid]->Fill(adc);}
00119 
00120     for (int j=96*2;j<96*3;j++){
00121       MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00122       zpos = channel->GetZ();
00123       hardid = channel->GetHardId();
00124       adc = channel->GetAnalogValue();
00125       hchamber_2[hardid]->Fill(adc);}
00126 
00127     for (int j=96*3;j<96*7;j++){
00128       MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00129       zpos = channel->GetZ();
00130       hardid = channel->GetHardId();
00131       adc = channel->GetAnalogValue();
00132       hchamber_3[hardid]->Fill(adc);}}}
00133 
00134 
00135 
00136 
00137 
00138 
00139 
00140 
00141   /*Fit*/
00142 
00143 
00144   cout<<endl;
00145   cout<<"Fitting"<<endl;
00146 
00147 
00148   TF1 * ff = new TF1("ff","gaus",0,1024);
00149   
00150   float param_chamber_0[96][3];
00151   float param_chamber_1[96][3];
00152   float param_chamber_2[96][3];
00153   float param_chamber_3[96*4][3];
00154 
00155 
00156   for (int i=0;i<96;i++){
00157 
00158     ff->SetParameters(hchamber_0[i]->GetMaximumBin(),hchamber_0[i]->GetMean(),hchamber_0[i]->GetRMS());
00159     hchamber_0[i]->Fit(ff,"q");
00160     param_chamber_0[i][0] = floor(0.5+ff->GetParameter(0));
00161     param_chamber_0[i][1] = floor(0.5+ff->GetParameter(1));
00162     param_chamber_0[i][2] = floor(0.5+ff->GetParameter(2));
00163 
00164     ff->SetParameters(hchamber_1[i]->GetMaximumBin(),hchamber_1[i]->GetMean(),hchamber_1[i]->GetRMS());
00165     hchamber_1[i]->Fit(ff,"q");
00166     param_chamber_1[i][0] = floor(0.5+ff->GetParameter(0));
00167     param_chamber_1[i][1] = floor(0.5+ff->GetParameter(1));
00168     param_chamber_1[i][2] = floor(0.5+ff->GetParameter(2));
00169 
00170     ff->SetParameters(hchamber_2[i]->GetMaximumBin(),hchamber_2[i]->GetMean(),hchamber_2[i]->GetRMS());
00171     hchamber_2[i]->Fit(ff,"q");
00172     param_chamber_2[i][0] = floor(0.5+ff->GetParameter(0));
00173     param_chamber_2[i][1] = floor(0.5+ff->GetParameter(1));
00174     param_chamber_2[i][2] = floor(0.5+ff->GetParameter(2));
00175 
00176     ff->SetParameters(hchamber_2[i]->GetMaximumBin(),hchamber_2[i]->GetMean(),hchamber_2[i]->GetRMS());
00177     hchamber_2[i]->Fit(ff,"q");
00178     param_chamber_2[i][0] = floor(0.5+ff->GetParameter(0));
00179     param_chamber_2[i][1] = floor(0.5+ff->GetParameter(1));
00180     param_chamber_2[i][2] = floor(0.5+ff->GetParameter(2));}
00181 
00182   for (int i=0;i<96*4;i++){
00183 
00184     ff->SetParameters(hchamber_3[i]->GetMaximumBin(),hchamber_3[i]->GetMean(),hchamber_3[i]->GetRMS());
00185     hchamber_3[i]->Fit(ff,"q");
00186     param_chamber_3[i][0] = floor(0.5+ff->GetParameter(0));
00187     param_chamber_3[i][1] = floor(0.5+ff->GetParameter(1));
00188     param_chamber_3[i][2] = floor(0.5+ff->GetParameter(2));}
00189 
00190 
00191 
00192 
00193   /*Writing*/
00194   cout<<endl;
00195   cout<<"Writing"<<endl;
00196   cout<<endl;
00197 
00198 //if (cut)
00199 //{ 
00200  ofstream out0(output_path+"PedestalADC1_cut.ini");
00201   ofstream out1(output_path+"PedestalADC2_cut.ini");
00202   ofstream out2(output_path+"PedestalADC3_cut.ini");
00203   ofstream out3(output_path+"PedestalADC4_cut.ini");
00204 
00205   out0<<"BLOCK0"<<endl;
00206   for (int i=0;i<96;i++){out0<<i<<","<<param_chamber_0[i][1]-pedestal<<","<<param_chamber_0[i][1]+param_chamber_0[i][2]*sigma_cut.Atoi()<<endl;}
00207   out0<<"-1"<<endl;
00208   out0<<"BLOCK1"<<endl;
00209   for (int i=0;i<96;i++){out0<<i<<","<<param_chamber_1[i][1]-pedestal<<","<<param_chamber_1[i][1]+param_chamber_1[i][2]*sigma_cut.Atoi()<<endl;}
00210   out0<<"-1"<<endl;
00211 
00212   out1<<"BLOCK0"<<endl;
00213   for (int i=0;i<96;i++){out1<<i<<","<<param_chamber_2[i][1]-pedestal<<","<<param_chamber_2[i][1]+param_chamber_2[i][2]*sigma_cut.Atoi()<<endl;}
00214   out1<<"-1"<<endl;
00215   out1<<"BLOCK1"<<endl;
00216   for (int i=0;i<96;i++){out1<<i<<","<<param_chamber_3[i][1]-pedestal<<","<<param_chamber_3[i][1]+param_chamber_3[i][2]*sigma_cut.Atoi()<<endl;}
00217   out1<<"-1"<<endl;
00218 
00219   out2<<"BLOCK0"<<endl;
00220   for (int i=0;i<96;i++){out2<<i<<","<<param_chamber_3[i+96][1]-pedestal<<","<<param_chamber_3[i][1]+param_chamber_3[i][2]*sigma_cut.Atoi()<<endl;}
00221   out2<<"-1"<<endl;
00222   out2<<"BLOCK1"<<endl;
00223   for (int i=0;i<96;i++){out2<<i<<","<<param_chamber_3[i+96*2][1]-pedestal<<","<<param_chamber_3[i][1]+param_chamber_3[i][2]*sigma_cut.Atoi()<<endl;}
00224   out2<<"-1"<<endl;
00225 
00226   out3<<"BLOCK0"<<endl;
00227   for (int i=0;i<96;i++){out3<<i<<","<<param_chamber_3[i+96*3][1]-pedestal<<","<<param_chamber_3[i][1]+param_chamber_3[i][2]*sigma_cut.Atoi()<<endl;}
00228   out3<<"-1"<<endl;
00229 //}
00230 //else
00231 //{
00232   ofstream out0(output_path+"PedestalADC1.ini");
00233   ofstream out1(output_path+"PedestalADC2.ini");
00234   ofstream out2(output_path+"PedestalADC3.ini");
00235   ofstream out3(output_path+"PedestalADC4.ini");
00236 
00237   out0<<"BLOCK0"<<endl;
00238   for (int i=0;i<96;i++){out0<<i<<","<<param_chamber_0[i][1]-pedestal<<","<<0<<endl;}
00239   out0<<"-1"<<endl;
00240   out0<<"BLOCK1"<<endl;
00241   for (int i=0;i<96;i++){out0<<i<<","<<param_chamber_1[i][1]-pedestal<<","<<0<<endl;}
00242   out0<<"-1"<<endl;
00243 
00244   out1<<"BLOCK0"<<endl;
00245   for (int i=0;i<96;i++){out1<<i<<","<<param_chamber_2[i][1]-pedestal<<","<<0<<endl;}
00246   out1<<"-1"<<endl;
00247   out1<<"BLOCK1"<<endl;
00248   for (int i=0;i<96;i++){out1<<i<<","<<param_chamber_3[i][1]-pedestal<<","<<0<<endl;}
00249   out1<<"-1"<<endl;
00250 
00251   out2<<"BLOCK0"<<endl;
00252   for (int i=0;i<96;i++){out2<<i<<","<<param_chamber_3[i+96][1]-pedestal<<","<<0<<endl;}
00253   out2<<"-1"<<endl;
00254   out2<<"BLOCK1"<<endl;
00255   for (int i=0;i<96;i++){out2<<i<<","<<param_chamber_3[i+96*2][1]-pedestal<<","<<0<<endl;}
00256   out2<<"-1"<<endl;
00257 
00258   out3<<"BLOCK0"<<endl;
00259   for (int i=0;i<96;i++){out3<<i<<","<<param_chamber_3[i+96*3][1]-pedestal<<","<<0<<endl;}
00260   out3<<"-1"<<endl;
00261 
00262 //}
00263   cout<<"Done."<<endl;
00264  
00265 }
00266
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/TB2010/Testbeam_june_10/Gassiplex/Pedestal_alignment/Pedestal_alignment.C
