Micromegas LAPP reconstruction

00001 #include <iostream>
00002 #include <map>
00003 #include "TKey.h"
00004 #include <TROOT.h>
00005 #include <TRint.h>
00006 #include <TKey.h>
00007 #include <TTree.h>
00008 #include <TFile.h>
00009 #include <TSystem.h>
00010 #include <TreeClass.hh>
00011 #include <TF1.h>
00012 #include <TH1I.h>
00013 #include <TH2I.h>
00014 #include <TH2F.h>
00015 #include <TGraph.h>
00016 #include <sstream>
00017 #include <vector>
00018 #include "Log.hh"
00019 
00020 #include "MicroException.hh"
00021 
00022 
00023 using namespace std;
00024 
00025 
00026 int main(int argc, char**argv){
00027 
00028   
00029   TString gassiplex_file = argv[1];
00030   TString hardroc2_file = argv[2];
00031   TString chip = argv[3];
00032   
00033   int tested_chip = chip.Atoi();
00034 
00035 
00036   TString dir = "/lapp_data/LC/Detecteurs/MicroMegas/data/TB2010/SPS_H4_june_2010/Root_files/";
00037   TString filename = dir + gassiplex_file + "_acq_HR2_" + hardroc2_file + "_merged.root";
00038 
00039   TString result_filename = "/LC/Detecteurs/MicroMegas/Offline/micromegasFrameWork/trunk/src/analyse/TB2010/Testbeam_june_10/Complete_detector/Results/" + gassiplex_file + "_acq_HR2_" + hardroc2_file + "_results.root";
00040   TFile * tf = new TFile(result_filename.Data(),"RECREATE");
00041 
00042 
00043 
00044   cout<<endl;
00045   cout<<"Gassiplex file "<<gassiplex_file<<endl;
00046   cout<<"Hardroc2 file  "<<hardroc2_file<<endl;
00047   cout<<"Tested chip    "<<chip<<endl;
00048   cout<<endl;
00049 
00050 
00051 
00052   bool verbose = 0;
00053 
00054 
00055   /*Variables*/
00056 
00057   vector <int> hot_event_index;
00058   vector <int> high_mult_event_index;
00059 
00060   TH2I * hxy[5];
00061   TH1I * hadc[4];
00062   TString name,title;
00063 
00064   for (int i=0;i<5;i++){
00065     name = Form("hxy_chb_%i",i+1);
00066     title = Form("Hit after cut in chamber %i;x (cm);y (cm)",i+1);
00067     hxy[i] = new TH2I(name,title,96,0,96,96,0,96);
00068     if (i!=4){
00069       name = Form("hadc_chb_%i",i+1);
00070       title = Form("ADC in chamber %i;ADC counts",i+1);
00071       hadc[i] = new TH1I(name,title,1024,0,1024);}}
00072 
00073   int nhigh_mult_max = 200;
00074   TH2I * hxy_m2_high_mult[nhigh_mult_max];
00075 
00076   for (int i=0;i<nhigh_mult_max;i++){
00077 
00078     name = Form("hxy_m2_high_mult_%i",i+1);
00079     title = Form("High multiplicity event number %i",i+1);
00080     hxy_m2_high_mult[i] = new TH2I(name,title,96,0,96,96,0,96);}
00081 
00082 
00083   TH1I * hmult_m2 = new TH1I("hmult_m2","Number of hits in m2",10,0,10);
00084   TH1I * hmult_m2_2 = new TH1I("hmult_m2_2","Number of hits in m2 for track extrapolated in a 2x2 cm2 region",10,0,10);
00085 
00086   TH1I * hntrack = new TH1I("hntrack","Number of telescope track extrapolated at the m2;region size (pad);number of track",10,0,10);
00087   TH1I * hntrack_m2 = new TH1I("hntrack_m2","Number of m2 track with extrapolation at the m2;region size (pad);number of track",10,0,10);
00088 
00089   TH2F * hxy_target = new TH2F("hxy_target","Position of telescope to m2 extrapolated hits;x (cm);y (cm);Ntest",96,0,96,96,0,96);
00090   TH1I * hdt = new TH1I("hdt","BCID distribution all chamber;BCID_dif - BCID_hit",400,0,400);
00091   TH1I * hdt_chip = new TH1I("hdt_chip","BCID distribution tested chip;BCID_dif - BCID_hit",400,0,400);
00092 
00093   TH1I * hdt_chip_channel[64];
00094   for (int i=0;i<64;i++){
00095     name = Form("hdt_chip_channel_%i",i);
00096     title = Form("BCID distribution tested chip of channel %i",i);
00097     hdt_chip_channel[i] = new TH1I(name,title,400,0,400);}
00098 
00099   TH1I * hdt_straight_tracks = new TH1I("hdt_chip","BCID distribution tested chip;BCID_dif - BCID_hit",400,0,400);
00100 
00101   TH2F * hxy_noise = new TH2F("hxy_noise","Expected number of noise hits in time window",8,0,8,8,0,8);
00102 
00103   TH1I * hnhit_m2_all_raw = new TH1I("hnhit_m2_all_raw","Number of hits per readout",1536,0,1536);
00104   TH1I * hnhit_m2_all_cut = new TH1I("hnhit_m2_all_cut","Number of hits per readout after cut",1536,0,1536);
00105   TH1I * hnhit_m2_chip_raw = new TH1I("hnhit_m2_chip_raw","Number of hits on chip per readout",64,0,64);
00106   TH1I * hnhit_m2_chip_cut = new TH1I("hnhit_m2_chip_cut","Number of hits on chip per readout after cut",64,0,64);
00107 
00108   TH2I * hxy_gassiplex = new TH2I("hxy_gassiplex","",96,0,96,96,0,96);
00109   TH2I * hxy_m2_all = new TH2I("hxy_m2_all","Hit in m2;x (cm);y (cm)",96,0,96,96,0,96);
00110   TH2I * hxy_m2_aligned = new TH2I("hxy_m2_aligned","Hit position in m2 for straight track;x (cm);y (cm)",96,0,96,96,0,96);
00111 
00112   TH2I * hxy_m2_aligned_max[3][3];
00113   
00114   for (int i=0;i<3;i++){
00115     for (int j=0;j<3;j++){
00116       name  = Form("hxy_m2_aligned_max_%i_%i",i,j);
00117       title = Form("Hit position in m2 for straight track going through profile max (%i,%i);x (cm);y (cm)",i-1,j-1);
00118       hxy_m2_aligned_max[i][j] = new TH2I(name,title,96,0,96,96,0,96);}}
00119 
00120   TH2I * hxy_m2_target = new TH2I("hxy_m2_target","Target hit position in m2;x (cm);y (cm)",96,0,96,96,0,96);
00121 
00122   TH1I * hnhot_chb = new TH1I("hnhot_chb","",6,0,6);
00123   TH1I * hnhot_chb_straight = new TH1I("hnhot_chb_straight","",6,0,6);
00124 
00125   TH1F * hx_residuals = new TH1F("hx_residuals","Residuals (extrapolated/m2 hit along x (2x2 region)",100,-10,10);
00126   TH1F * hy_residuals = new TH1F("hy_residuals","Residuals (extrapolated/m2 hit along y (2x2 region)",100,-10,10);
00127 
00128   TH1F * hx_m2 = new TH1F("hx_m2","",96,0,96);
00129   TH1F * hy_m2 = new TH1F("hy_m2","",96,0,96);
00130   TH1F * hx_m2_all = new TH1F("hx_m2_all","",96,0,96);
00131   TH1F * hy_m2_all = new TH1F("hy_m2_all","",96,0,96);
00132 
00133   TH1F * hpres = new TH1F("hpres","Pressure",500,900,1000);
00134   TH1F * htemp = new TH1F("htemp","Temperature",500,280,310);
00135 
00136 
00137   /*Cuts*/
00138 
00139   int adc_cut = 35;
00140   int dt_cut_min = 6;
00141   int dt_cut_max = 7;
00142 
00143   Long64_t dt = 0;
00144 
00145   int threshold = 0;
00146   int nchannel = 0;
00147   int chamber_id = 0;
00148   int hard_id = 0;
00149   int chip_id = 0;
00150   int channel_hard_id = 0;
00151   int adc = 0;
00152   int xpos = 0;
00153   int ypos = 0;
00154   int nhot_gassiplex = 0;
00155   int nhot_m2 = 0;
00156   int nhot_m2_2 = 0;
00157   int nhot_m2_4 = 0;
00158   int nhot_m2_6 = 0;
00159   int nhot_m2_8 = 0;
00160 
00161   bool channel_set[64];
00162   int channel_x[64];
00163   int channel_y[64];
00164   for (int i=0;i<64;i++){
00165     channel_x[i] = 0;
00166     channel_y[i] = 0;
00167     channel_set[i] = 0;}
00168 
00169   int nhit_chb[5];
00170   int nhit_m2_raw = 0;
00171   int nhit_m2_chip_raw = 0;
00172   int nhit_m2_chip_cut = 0;
00173 
00174   int nhot_event = 0;
00175   int nhot_m2_chip = 0;
00176   int ntest = 0;
00177   int nhit = 0;
00178   int ntest_2 = 0;
00179   int ntest_4 = 0;
00180   int ntest_6 = 0;
00181   int ntest_8 = 0;
00182   int nout_chip = 0;
00183   int intime_hit = 0;
00184   int noise_hit = 0;
00185   int noise_hit_2 = 0;
00186   int noise_hit_4 = 0;
00187   int noise_hit_6 = 0;
00188   int noise_hit_8 = 0;
00189   int nmult = 0;
00190   int nmult_2 = 0;
00191 
00192   double signal_to_noise = 0;
00193   double noise_density = 0;
00194   double noise_probability = 0;
00195   double eff_2 = 0;
00196   double eff_4 = 0;
00197   double eff_6 = 0;
00198   double eff_8 = 0;
00199 
00200   bool hot_chb[5];
00201   bool hit = 0;
00202   bool straight_track = 0;
00203   bool hot_event = 0;
00204   bool target_m2_2_bis = 0;
00205   bool target_m2_2 = 0;
00206   bool target_m2_4 = 0;
00207   bool target_m2_6 = 0;
00208   bool target_m2_8 = 0;
00209   bool hit_already = 0;
00210 
00211   int xtarget_int = 0;
00212   int ytarget_int = 0;
00213   int delta_max_x = 0;
00214   int delta_max_y = 0;
00215   double xtarget = 0;
00216   double ytarget = 0;
00217   double xpos_telescope = 0;
00218   double ypos_telescope = 0;
00219 
00220   TF1 * fgaus_x = new TF1("fgaus_x","gaus",0,96);
00221   TF1 * fgaus_y = new TF1("fgaus_y","gaus",0,96);
00222 
00223   double par0[5] = {0,0,0,0,0};
00224   double par1_x[5] = {0,0,0,0,0};
00225   double par1_y[5] = {0,0,0,0,0};
00226   double par2_x[5] = {0,0,0,0,0};
00227   double par2_y[5] = {0,0,0,0,0};
00228   double par1_err_x[5] = {0,0,0,0,0};
00229   double par1_err_y[5] = {0,0,0,0,0};
00230   double par2_err_x[5] = {0,0,0,0,0};
00231   double par2_err_y[5] = {0,0,0,0,0};
00232   double chi2_ndf_x[5] = {0,0,0,0,0};
00233   double chi2_ndf_y[5] = {0,0,0,0,0};
00234 
00235   double integral = 0;
00236   double mean_x = 0;
00237   double mean_y = 0;
00238   double sigma_x = 0;
00239   double sigma_y = 0;
00240 
00241   double gassiplex_mean_x = 0;
00242   double gassiplex_mean_y = 0;
00243   double gassiplex_mean_err_x = 0;
00244   double gassiplex_mean_err_y = 0;
00245 
00246   double gassiplex_sigma_x = 0;
00247   double gassiplex_sigma_y = 0;
00248   double gassiplex_sigma_err_x = 0;
00249   double gassiplex_sigma_err_y = 0;
00250 
00251   double gassiplex_profile_max_x = 0;
00252   double gassiplex_profile_max_y = 0;
00253 
00254   double m2_profile_max_x = 0;
00255   double m2_profile_max_y = 0;
00256 
00257   double m2_mean_x = 0;
00258   double m2_mean_y = 0;
00259   double m2_mean_err_x = 0;
00260   double m2_mean_err_y = 0;
00261 
00262   double m2_sigma_x = 0;
00263   double m2_sigma_y = 0;
00264   double m2_sigma_err_x = 0;
00265   double m2_sigma_err_y = 0;
00266 
00267   double offset_x = 0;
00268   double offset_y = 0;
00269   double offset_err_x = 0;
00270   double offset_err_y = 0;
00271 
00272   double xmin_chip = 0;
00273   double ymin_chip = 0;
00274   double xmax_chip = 0;
00275   double ymax_chip = 0;
00276 
00277 
00278 
00279   TFile f(filename);
00280   cout<<endl;cout<<"reading File : "<<filename<<endl;cout<<endl;
00281 
00282   TIter nextkey(f.GetListOfKeys());
00283   TKey *key;
00284 
00285   int nevent = 0;
00286 
00287   TTree * tree;
00288   MTRun * run;
00289   MTEvent *evt = new MTEvent();
00290 
00291 
00292 
00293 
00294   while (key = (TKey*)nextkey()) {
00295 
00296     tree = (TTree*)key->ReadObj();                
00297     run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
00298     run->Info();
00299 
00300 
00301 
00302     /*Get chip thresholds*/
00303 
00304     for( map<Long_t,MTChip*>::const_iterator iiChip=run->GetChipConfigs().begin(); iiChip!=run->GetChipConfigs().end(); ++iiChip){
00305       const MTChip& chip  = *((*iiChip).second);
00306       try{
00307         const MTHardroc2Chip &hr2Chip = dynamic_cast<const MTHardroc2Chip &> (chip);
00308         threshold = hr2Chip.GetThresholdDac_0();
00309         cout <<"Chip "<<chip.GetId()<<" threshold = "<<hr2Chip.GetThresholdDac_0()<<" "<<hr2Chip.GetThresholdDac_1()<<" "<<hr2Chip.GetThresholdDac_2()<<endl;}
00310       catch (...) {}}
00311 
00312 
00313 
00314 
00315     
00316     nevent = tree->GetEntries();
00317     cout <<"Total number of event (or trigger) ="<<nevent<<endl;
00318     
00319     TBranch *branch= tree->GetBranch("MTEvent");
00320     branch->SetAddress(&evt);
00321 
00322 
00323 
00324 
00325 
00326   /*Determine beam profile in each Gassiplex chamber and in m2*/
00327   /*Store index of events with straight tracks (single aligned hits in gassiplex)*/
00328 
00329     for (int i=0;i<nevent;i++){
00330     //for (int i=0;i<500;i++){
00331 
00332     hxy_gassiplex->Reset();
00333     nhot_m2 = 0;
00334     nhot_gassiplex = 0;
00335     nhit_m2_chip_raw = 0;
00336     nhit_m2_chip_cut = 0;
00337     nhit_m2_raw = 0;
00338     for (int k=0;k<5;k++){hot_chb[k] = 0;nhit_chb[k] = 0;}
00339 
00340     tree->GetEntry(i);
00341     nchannel = evt->GetNchannel();    
00342 
00343     hpres->Fill(evt->GetPressure());
00344     htemp->Fill(evt->GetTemperature()+273);
00345 
00346 
00347     for (int j=0;j<nchannel;j++){
00348       
00349       hit = 0;
00350 
00351       MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00352       
00353       chamber_id = channel->GetChamberId();
00354 
00355       /*Gassiplex*/
00356 
00357       if (chamber_id>=1 && chamber_id<=4){
00358         adc = channel->GetAnalogValue();
00359         channel_hard_id = channel->GetHardId();
00360         hadc[chamber_id-1]->Fill(adc);
00361         if (adc>=adc_cut){
00362           hit = 1;
00363           hot_chb[chamber_id-1] = 1;
00364           xpos = channel->GetX();
00365           ypos = channel->GetY();
00366           if (verbose){cout<<"gassiplex hit at "<<xpos<<" "<<ypos<<" in chamber "<<chamber_id<<endl;}
00367           hxy_gassiplex->Fill(xpos,ypos);
00368           hxy[chamber_id-1]->Fill(xpos,ypos);
00369           if (nhit_chb[chamber_id-1]==0){nhot_gassiplex++;}
00370           nhit_chb[chamber_id-1]++;}}
00371 
00372       /*M2 chamber*/
00373 
00374       if (chamber_id==5){       
00375 
00376         nhit_m2_raw++;
00377         chip_id = channel->GetChipId();
00378         if (chip_id==tested_chip){nhit_m2_chip_raw++;}
00379 
00380         dt = channel->GetBcId_Dif() - channel->GetBcId_Hit();
00381         xpos = channel->GetX();
00382         ypos = channel->GetY();
00383 
00384         hdt->Fill(dt);
00385         hx_m2_all->Fill(xpos);
00386         hy_m2_all->Fill(ypos);
00387         hxy_m2_all->Fill(xpos,ypos);
00388 
00389         if (channel->GetChipId()==tested_chip){
00390 
00391           hard_id = channel->GetHardId();
00392 
00393           if (!channel_set[hard_id]){       
00394             channel_x[hard_id]=xpos%8;
00395             channel_y[hard_id]=ypos%8;
00396             channel_set[hard_id] = 1;}
00397 
00398           hdt_chip_channel[hard_id]->Fill(dt);
00399           hdt_chip->Fill(dt);}
00400 
00401         if (dt>=dt_cut_min && dt<=dt_cut_max){
00402 
00403           hit = 1;
00404           hot_chb[chamber_id-1] = 1;
00405           if (verbose){cout<<"m2 chamber hit at "<<xpos<<" "<<ypos<<endl;}
00406           hxy[chamber_id-1]->Fill(xpos,ypos);
00407           hx_m2->Fill(xpos);
00408           hy_m2->Fill(ypos);
00409           nhit_chb[chamber_id-1]++;
00410           if (chip_id==tested_chip){nhit_m2_chip_cut++;}
00411           nhot_m2 = 1;}}}
00412 
00413     hnhit_m2_all_raw->Fill(nhit_m2_raw);
00414     hnhit_m2_all_cut->Fill(nhit_chb[4]);
00415     hnhit_m2_chip_raw->Fill(nhit_m2_chip_raw);
00416     hnhit_m2_chip_cut->Fill(nhit_m2_chip_cut);
00417 
00418     /*Straight tracks if single aligned hits in at least 3 chambers*/
00419 
00420     straight_track = (hxy_gassiplex->GetRMS(1)==0) && (hxy_gassiplex->GetRMS(2)==0);
00421     hot_event = (nhot_gassiplex==3) || (nhot_gassiplex==4);
00422 
00423     if (verbose){
00424       cout<<"Event "<<i+1<<" -> ";
00425       for (int k=0;k<5;k++){cout<<nhit_chb[k]<<" ";}cout<<endl;}
00426 
00427     if (straight_track){hnhot_chb_straight->Fill(nhot_gassiplex+nhot_m2);}
00428     if (hot_event && straight_track){hot_event_index.push_back(i);}}}
00429 
00430 
00431 
00432 
00433 
00434 
00435 
00436   /*Calculate average pressure and temperature*/
00437 
00438   cout<<endl;
00439   cout<<"Pressure = "<<hpres->GetMean()<<" +/- "<<hpres->GetRMS()<<" mbar"<<endl;
00440   cout<<"Temperature = "<<htemp->GetMean()<<" +/- "<<htemp->GetRMS()<<" K"<<endl;
00441   cout<<endl;
00442 
00443 
00444 
00445 
00446   /*Fit beam profiles*/
00447 
00448   for (int i=0;i<5;i++){
00449 
00450     integral = hxy[i]->GetEntries();
00451     mean_x   = hxy[i]->GetMean(1);
00452     mean_y   = hxy[i]->GetMean(2);
00453     sigma_x  = hxy[i]->GetRMS(1);
00454     sigma_y  = hxy[i]->GetRMS(2);
00455 
00456     fgaus_x->SetParameters(integral,mean_x,sigma_x);
00457     fgaus_y->SetParameters(integral,mean_y,sigma_y);
00458 
00459     hxy[i]->ProjectionX()->Fit(fgaus_x,"q");
00460     hxy[i]->ProjectionY()->Fit(fgaus_y,"q");
00461 
00462     par0[i] = integral;
00463     par1_x[i] = fgaus_x->GetParameter(1);
00464     par1_y[i] = fgaus_y->GetParameter(1);
00465     par2_x[i] = fgaus_x->GetParameter(2);
00466     par2_y[i] = fgaus_y->GetParameter(2);
00467 
00468     par1_err_x[i] = fgaus_x->GetParError(1);
00469     par1_err_y[i] = fgaus_y->GetParError(1);
00470     par2_err_x[i] = fgaus_x->GetParError(2);
00471     par2_err_y[i] = fgaus_y->GetParError(2);
00472 
00473     chi2_ndf_x[i] = fgaus_x->GetChisquare()/fgaus_x->GetNDF();
00474     chi2_ndf_y[i] = fgaus_y->GetChisquare()/fgaus_y->GetNDF();
00475     
00476     cout<<endl;
00477     cout<<"Chamber "<<i+1<<" -> "<<par0[i]<<" hits"<<endl;
00478     cout<<" mean     = ("<<par1_x[i]<<","<<par1_y[i]<<") +/- ("<<par1_err_x[i]<<","<<par1_err_y[i]<<")"<<endl;
00479     cout<<" sigma    = ("<<par2_x[i]<<","<<par2_y[i]<<") +/- ("<<par2_err_x[i]<<","<<par2_err_y[i]<<")"<<endl;
00480     cout<<" chi2/ndf = ("<<chi2_ndf_x[i]<<","<<chi2_ndf_y[i]<<")"<<endl;
00481 
00482     if (i<=3){
00483       
00484       gassiplex_mean_x += par1_x[i]/4.;
00485       gassiplex_mean_y += par1_y[i]/4.;
00486 
00487       gassiplex_mean_err_x += par1_err_x[i]/4.;
00488       gassiplex_mean_err_y += par1_err_y[i]/4.;
00489 
00490       gassiplex_sigma_x += par2_x[i]/4.;
00491       gassiplex_sigma_y += par2_y[i]/4.;
00492 
00493       gassiplex_sigma_err_x += par2_err_x[i]/4.;
00494       gassiplex_sigma_err_y += par2_err_y[i]/4.;}}
00495 
00496   gassiplex_profile_max_x = floor(gassiplex_mean_x);
00497   gassiplex_profile_max_y = floor(gassiplex_mean_y);
00498 
00499 
00500   m2_mean_x = par1_x[4];
00501   m2_mean_y = par1_y[4];
00502 
00503   m2_mean_err_x = par1_err_x[4];
00504   m2_mean_err_y = par1_err_y[4];
00505 
00506   m2_sigma_x = par2_x[4];
00507   m2_sigma_y = par2_y[4];
00508 
00509   m2_sigma_err_x = par2_err_x[4];
00510   m2_sigma_err_y = par2_err_y[4];
00511 
00512   m2_profile_max_x = floor(m2_mean_x);
00513   m2_profile_max_y = floor(m2_mean_y);  
00514 
00515   offset_x = gassiplex_mean_x - m2_mean_x;
00516   offset_y = gassiplex_mean_y - m2_mean_y;
00517   offset_err_x = gassiplex_mean_err_x + m2_mean_err_x;
00518   offset_err_y = gassiplex_mean_err_y + m2_mean_err_y;
00519   
00520   cout<<endl;cout<<"Beam profile"<<endl;cout<<endl;
00521   cout<<" Gassiplex telescope"<<endl;
00522   cout<<"   ("<<gassiplex_mean_x<<","<<gassiplex_mean_y<<") +/- ("<<gassiplex_mean_err_x<<","<<gassiplex_mean_err_y<<")"<<endl;
00523   cout<<" M2 telescope"<<endl;
00524   cout<<"   ("<<m2_mean_x<<","<<m2_mean_y<<") +/- ("<<m2_mean_err_x<<","<<m2_mean_err_y<<")"<<endl;
00525   cout<<" Offset (x,y)_telescope -> -Delta (x,y) -> (x,y)_m2"<<endl;
00526   cout<<"   ("<<offset_x<<","<<offset_y<<") +/- ("<<offset_err_x<<","<<offset_err_y<<")"<<endl;
00527   cout<<endl;
00528 
00529 
00530 
00531   cout<<"Tested chip number is "<<tested_chip<<endl;
00532 
00533   xmin_chip = floor(m2_profile_max_x/8.) * 8;
00534   xmax_chip = xmin_chip + 8;
00535 
00536   ymin_chip = floor(m2_profile_max_y/8.) * 8;
00537   ymax_chip = ymin_chip + 8;
00538 
00539   cout<<" xmin = "<<xmin_chip<<"  xmax = "<<xmax_chip<<endl;
00540   cout<<" ymin = "<<ymin_chip<<"  ymax = "<<ymax_chip<<endl;cout<<endl;
00541 
00542 
00543 
00544 
00545 
00546 
00547 
00548   /*Loop over straight track events*/
00549 
00550 
00551   nhot_m2 = 0;
00552   nhot_event = hot_event_index.size();
00553   cout<<nevent<<" events in total"<<endl;
00554   cout<<"Loop over "<<nhot_event<<" hot events (single aligned hits in at least 3 chambers)"<<endl;
00555 
00556 
00557   
00558   for (int i=0;i<nhot_event;i++){  
00559     
00560     tree->GetEntry(hot_event_index[i]);
00561     nchannel = evt->GetNchannel();
00562     
00563     
00564     /*Get coordinates of the gassiplex hit*/
00565     /*transform them into m2 frame*/
00566     
00567     for (int j=0;j<nchannel;j++){
00568       
00569       MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00570       chamber_id = channel->GetChamberId();
00571 
00572       
00573       /*Extrapolate gassiplex hits to the m2 chamber*/
00574 
00575       if (chamber_id==1 || chamber_id==2 || chamber_id==3 || chamber_id==4){
00576         adc = channel->GetAnalogValue();
00577         if (adc>adc_cut){
00578 
00579           xpos_telescope = channel->GetX();
00580           ypos_telescope = channel->GetY();
00581 
00582           xtarget = xpos_telescope - offset_x;
00583           ytarget = ypos_telescope - offset_y;
00584 
00585           xtarget_int = floor(xtarget);
00586           ytarget_int = floor(ytarget);
00587 
00588           hxy_target->Fill(xtarget,ytarget);
00589 
00590           /*2x2 region*/
00591           if (((xtarget_int>=(xmin_chip+3))&&(xtarget_int<=(xmin_chip+4)) &&
00592                (ytarget_int>=(ymin_chip+3))&&(ytarget_int<=(ymin_chip+4)))){
00593             target_m2_2 = 1;target_m2_4 = 1;target_m2_6 = 1;target_m2_8 = 1;
00594             ntest_2++;ntest_4++;ntest_6++;ntest_8++;}
00595 
00596           else{
00597 
00598             /*4x4 region*/
00599             if (((xtarget_int>=(xmin_chip+2))&&(xtarget_int<=(xmin_chip+5)) &&
00600                  (ytarget_int>=(ymin_chip+2))&&(ytarget_int<=(ymin_chip+5)))){
00601               target_m2_4 = 1;target_m2_6 = 1;target_m2_8 = 1;
00602               ntest_4++;ntest_6++;ntest_8++;}
00603 
00604             else{
00605 
00606               /*6x6 region*/
00607               if (((xtarget_int>=(xmin_chip+1))&&(xtarget_int<=(xmin_chip+6)) &&
00608                    (ytarget_int>=(ymin_chip+1))&&(ytarget_int<=(ymin_chip+6)))){
00609                 target_m2_6 = 1;target_m2_8 = 1;
00610                 ntest_6++;ntest_8++;}
00611 
00612               else{
00613 
00614                 /*8x8 region*/
00615                 if (((xtarget_int>=xmin_chip)&&(xtarget_int<=xmin_chip+7) &&
00616                      (ytarget_int>=ymin_chip)&&(ytarget_int<=ymin_chip+7))){
00617                   target_m2_8 = 1;
00618                   ntest_8++;}
00619 
00620                 else{nout_chip++;}}}}
00621 
00622           j = nchannel;}}}
00623 
00624 
00625     
00626     /*Look for hits in m2 chamber*/
00627     
00628     target_m2_2_bis = target_m2_2;
00629 
00630       for (int j=0;j<nchannel;j++){
00631         
00632         MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00633         chamber_id = channel->GetChamberId();
00634         
00635         if (chamber_id==5){     
00636           dt = channel->GetBcId_Dif() - channel->GetBcId_Hit();
00637           if (dt>=dt_cut_min && dt<=dt_cut_max){
00638 
00639             nhot_m2++;
00640             xpos = channel->GetX();
00641             ypos = channel->GetY();
00642             hxy_m2_target->Fill(xtarget,ytarget);
00643             hxy_m2_aligned->Fill(xpos,ypos);
00644 
00645             delta_max_x = xpos_telescope-gassiplex_profile_max_x;
00646             delta_max_y = ypos_telescope-gassiplex_profile_max_y;
00647 
00648             if ((fabs(delta_max_x)<=1)&&(fabs(delta_max_y)<=1)){hxy_m2_aligned_max[delta_max_x+1][delta_max_y+1]->Fill(xpos,ypos);}
00649 
00650             if (channel->GetChipId() == tested_chip){
00651 
00652               nmult++;
00653 
00654               if (target_m2_2_bis){
00655 
00656                 nmult_2++;
00657                 hx_residuals->Fill(xtarget - xpos);
00658                 hy_residuals->Fill(ytarget - ypos);}
00659 
00660               if (hit_already==0){
00661 
00662                 nhot_m2_chip++;
00663                 
00664                 if (target_m2_2){nhot_m2_2++;}
00665                 if (target_m2_4){nhot_m2_4++;}
00666                 if (target_m2_6){nhot_m2_6++;}
00667                 if (target_m2_8){nhot_m2_8++;}
00668 
00669                 hit_already = 1;}}}}}
00670 
00671       if (nmult>1){high_mult_event_index.push_back(hot_event_index[i]);}
00672       if (nmult==5){cout<<"event "<<hot_event_index[i]<<" has 5 hits in time"<<endl;}
00673 
00674       target_m2_2 = 0;
00675       target_m2_4 = 0;
00676       target_m2_6 = 0;
00677       target_m2_8 = 0;
00678 
00679       hmult_m2->Fill(nmult);
00680       hmult_m2_2->Fill(nmult_2);
00681       nmult = 0;
00682       nmult_2 = 0;
00683       hit_already = 0;}
00684   
00685 
00686 
00687   hntrack->SetBinContent(2+1,nhot_m2_2);
00688   hntrack->SetBinContent(4+1,nhot_m2_4);
00689   hntrack->SetBinContent(6+1,nhot_m2_6);
00690   hntrack->SetBinContent(8+1,nhot_m2_8);
00691 
00692   hntrack_m2->SetBinContent(2+1,ntest_2);
00693   hntrack_m2->SetBinContent(4+1,ntest_4);
00694   hntrack_m2->SetBinContent(6+1,ntest_6);
00695   hntrack_m2->SetBinContent(8+1,ntest_8);
00696 
00697 
00698 
00699 
00700 
00701 
00702 
00703 
00704   /*Save events with high multiplicity (>1)*/
00705 
00706   int nhigh_mult_event = high_mult_event_index.size();
00707 
00708   cout<<endl;
00709   cout<<nhigh_mult_event<<" found with multiplicity larger than 1"<<endl;
00710   cout<<endl;
00711   
00712 
00713   /*Loop over events*/
00714 
00715   for (int i=0;i<nhigh_mult_event;i++){  
00716 
00717     if (i<nhigh_mult_max){
00718     
00719     tree->GetEntry(high_mult_event_index[i]);
00720     nchannel = evt->GetNchannel();
00721     
00722     for (int j=0;j<nchannel;j++){
00723       
00724       MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00725       chamber_id = channel->GetChamberId();
00726 
00727       if (chamber_id==5){
00728 
00729         dt = channel->GetBcId_Dif() - channel->GetBcId_Hit();
00730           if (dt>=dt_cut_min && dt<=dt_cut_max){
00731 
00732             xpos = channel->GetX();
00733             ypos = channel->GetY();
00734 
00735             title = Form("High multiplicity event number %i",high_mult_event_index[i]);
00736             hxy_m2_high_mult[i]->SetTitle(title);
00737 
00738             hxy_m2_high_mult[i]->Fill(xpos,ypos);}}}}}
00739 
00740 
00741 
00742 
00743   /*Noise hit calculation*/
00744   
00745   cout<<endl;
00746   cout<<"Noise per channel of tested chip"<<endl;
00747   cout<<endl;
00748 
00749   float noise_sum = 0;
00750 
00751   for (int i=0;i<64;i++){
00752 
00753     noise_density = hdt_chip_channel[i]->Integral(10,400)/390.;
00754     cout<<"Channel "<<i<<" at "<<channel_x[i]<<" "<<channel_y[i]<<" -> "<<noise_density * (dt_cut_max - dt_cut_min+1)<<endl;
00755     noise_sum += noise_density * (dt_cut_max - dt_cut_min+1);
00756     hxy_noise->SetBinContent(channel_x[i]+1,channel_y[i]+1,noise_density * (dt_cut_max - dt_cut_min+1));}
00757 
00758   cout<<" in total -> "<<noise_sum<<" hits"<<endl;
00759   cout<<endl;
00760 
00761 
00762 
00763   noise_density = hdt_chip->Integral(10,400)/390.;
00764   noise_hit = noise_density * (dt_cut_max - dt_cut_min+1);
00765   intime_hit = hdt_chip->Integral(dt_cut_min+1,dt_cut_max+1);
00766   signal_to_noise = intime_hit/(1.*noise_hit);
00767   noise_probability = noise_hit/(1.*intime_hit);
00768 
00769   noise_hit_2 = floor(nhot_m2_2 * noise_probability) + 1;
00770   noise_hit_4 = floor(nhot_m2_4 * noise_probability) + 1;
00771   noise_hit_6 = floor(nhot_m2_6 * noise_probability) + 1;
00772   noise_hit_8 = floor(nhot_m2_8 * noise_probability) + 1;
00773 
00774   eff_2 = (nhot_m2_2 - noise_hit_2)/(1.*ntest_2)*100.;
00775   eff_4 = (nhot_m2_4 - noise_hit_4)/(1.*ntest_4)*100.;
00776   eff_6 = (nhot_m2_6 - noise_hit_6)/(1.*ntest_6)*100.;
00777   eff_8 = (nhot_m2_8 - noise_hit_8)/(1.*ntest_8)*100.;
00778 
00779 
00780   cout<<endl;
00781   cout<<nhot_event<<" tracks found in telescope"<<endl;
00782   cout<<"  "<<ntest_2<<" extrapolate to a 2x2 region"<<endl;
00783   cout<<"  "<<ntest_4<<" extrapolate to a 4x4 region"<<endl;
00784   cout<<"  "<<ntest_6<<" extrapolate to a 6x6 region"<<endl;
00785   cout<<"  "<<ntest_8<<" extrapolate to a 8x8 region"<<endl;
00786   cout<<"  "<<nout_chip<<" extrapolate outside tested chip"<<endl;
00787   cout<<endl;
00788   cout<<nhot_m2<<" tracks found in m2"<<endl;
00789   cout<<"  "<<nhot_m2_chip<<" track found in tested chip = "<<nhot_m2_chip/(1.*nhot_m2)*100.<<" %"<<endl;
00790   cout<<endl;
00791 
00792   cout<<"Signal to noise ratio = ("<<intime_hit<<" - "<<noise_hit<<") / "<<noise_hit<<" = "<<intime_hit-noise_hit<<" / "<<noise_hit<<" = "<<(intime_hit-noise_hit)/(1.*noise_hit)<<endl;
00793   cout<<"Noise hit probability = "<<noise_hit<<" / "<<intime_hit<<" = "<<noise_probability * 100.<<" %"<<endl;
00794   cout<<endl;
00795 
00796   cout<<"Region \t\t NTrack / (Nhit_m2 - Nhit_noise) = \t Efficiency"<<endl;
00797 
00798   cout<<"2 x 2 cm2 \t "<<ntest_2<<" / ("<<nhot_m2_2<<" - "<<noise_hit_2<<") = \t\t\t "<<eff_2<<" %"<<endl;
00799   cout<<"4 x 4 cm2 \t "<<ntest_4<<" / ("<<nhot_m2_4<<" - "<<noise_hit_4<<") = \t\t\t "<<eff_4<<" %"<<endl;
00800   cout<<"6 x 6 cm2 \t "<<ntest_6<<" / ("<<nhot_m2_6<<" - "<<noise_hit_6<<") = \t\t\t "<<eff_6<<" %"<<endl;
00801   cout<<"8 x 8 cm2 \t "<<ntest_8<<" / ("<<nhot_m2_8<<" - "<<noise_hit_8<<") = \t\t\t "<<eff_8<<" %"<<endl;
00802   
00803 
00804 
00805 
00806 
00807   tf->Write();
00808   tf->Close();
00809 
00810 
00811 
00812 
00813 
00814 
00815 
00816 return 0;}
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/TB2010/Testbeam_june_10/Complete_detector/Analysis_tb_m2_2.C
