Micromegas LAPP reconstruction

00001  {
00002    gROOT->Reset();
00003    gSystem->Load("/LC/Detecteurs/MicroMegas/Offline/micromegasFrameWork/trunk/lib/libMicro.so");
00004 
00005 
00006 
00007    int adc_cut = 30; 
00008 
00009    bool fit = 0;
00010    bool print = 0;
00011 
00012    bool rotate90 = 1;
00013 
00014 
00015    //READ INFO FROM XML
00016    std::ifstream xml;
00017    xml.open("telescope_beam_profile.xml",ifstream::in);
00018    std::string line = "";
00019    std::string lineb = "";
00020    while(line.find("path=")==string::npos || lineb.find("<output")==string::npos)
00021      {
00022         lineb=line;
00023         xml>>line;
00024      }
00025    int start = line.find("\"")+1;
00026    int end =  line.rfind("\"")-6;
00027    TString file = line.substr(start,end);
00028    //TString file = "p2507_2.root";    
00029 
00030 
00031    file = "/gpfs/LAPP-DATA/LC/Detecteurs/MicroMegas/data/TB2011/SPS_H4_aug_2011/Telescope/Beam/Root_files/b0708_4.root";
00032 
00033 
00034 
00035 
00036    bool noisy[9][96];
00037    for (int i=0;i<9;i++){
00038      for (int j=0;j<96;j++){
00039          noisy[i][j] = 0;}}
00040 
00041    /*HV pad*/
00042    noisy[6][77] = 1;
00043    noisy[7][77] = 1;
00044    noisy[8][77] = 1;
00045    
00046    /*floating pad*/
00047    noisy[6][9] = 1;
00048    noisy[8][33] = 1;
00049 
00050 
00051 
00052 
00053 
00054 
00055 
00056    TH1I * hc_adc[9];
00057    TH1I * hs_adc[9];
00058    TH2I * hxy_hit[9];
00059    TH2I * hxy_adc[9];
00060    TH1I * hadc[9][96];
00061    TH1I * hnchan[9];
00062    TH1I * hnhit[9];
00063    TH1I * hmult[9];
00064    TF1 *  fadc[9][96];
00065    TF1 *  fx[9];
00066    TF1 *  fy[9];
00067 
00068    int nhit[9];
00069    int nchan[9];
00070 
00071    TH2I * hxy_pad_temp = new TH2I("hxy_pad_temp","",32,0,32,32,0,32);
00072    TH1I * hnhot = new TH1I("hnhot","Number of hot telescope chambers (Nhit>0);Nhot",9,0,9);
00073    TH1I * hnsuper_hot = new TH1I("hnsuper_hot","Number of super hot telescope chambers (Nhit=1);Nhot",9,0,9);
00074 
00075 
00076    TString name,title;
00077         
00078    for (int i=0;i<9;i++)
00079      {
00080 
00081        nhit[i] = 0;
00082        nchan[i] = 0;
00083 
00084        name = Form("hc_adc_chamber_%i",i+1);
00085        title = Form("ADC pattern chamber %i;channel;ADC counts",i+1);
00086        hc_adc[i] = new TH1I(name,title,96,0,96);
00087 
00088        name = Form("hs_adc_chamber_%i",i+1);
00089        title = Form("ADC pattern chamber %i;strip number;ADC counts",i+1);
00090        hs_adc[i] = new TH1I(name,title,96,0,96);
00091 
00092        name = Form("hxy_hit_chamber_%i",i+1);
00093        title = Form("Hit pattern (adc>%i) chamber %i;y (cm);x (cm)",adc_cut,i+1);
00094        hxy_hit[i] = new TH2I(name,title,32,0,32,32,0,32);
00095 
00096        name = Form("hxy_adc_chamber_%i",i+1);
00097        title = Form("ADC sum pattern chamber %i;y (cm);x (cm)",adc_cut,i+1);
00098        hxy_adc[i] = new TH2I(name,title,32,0,32,32,0,32);
00099 
00100        name = Form("hnchan_chamber_%i",i+1);
00101        title = Form("Number of readout channel - chamber %i;channel number",i+1);
00102        hnchan[i] = new TH1I(name,title,97,0,97);
00103        hnchan[i]->SetFillColor(2);
00104 
00105        name = Form("hnhit_chamber_%i",i+1);
00106        title = Form("Number of hit (adc>%i) - chamber %i;channel number",adc_cut,i+1);
00107        hnhit[i] = new TH1I(name,title,97,0,97);
00108        hnhit[i]->SetFillColor(3);
00109 
00110        name = Form("hmult_chamber_%i",i+1);
00111        title = Form("Hit multiplicity (adc>%i) - chamber %i;channel number",adc_cut,i+1);
00112        hmult[i] = new TH1I(name,title,97,0,97);
00113        hmult[i]->SetFillColor(3);
00114 
00115        name = Form("fx_chamber_%i",i+1);
00116        fx[i] = new TF1(name,"gaus",0,32);
00117        fx[i]->SetLineColor(1);
00118        fx[i]->SetLineWidth(2);
00119 
00120        name = Form("fy_chamber_%i",i+1);
00121        fy[i] = new TF1(name,"gaus",0,32);
00122        fy[i]->SetLineColor(1);
00123        fy[i]->SetLineWidth(2);
00124 
00125        for (int j=0;j<96;j++)
00126          {
00127            name = Form("hadc_chamber_%i_channel_%i",i+1,j);
00128            title = Form("ADC of chamber %i channel %i",i+1,j);
00129            hadc[i][j] = new TH1I(name,title,1024,0,1024);
00130            hadc[i][j]->SetFillColor(2);
00131 
00132            name = Form("fadc_chamber_%i_channel_%i",i+1,j);
00133            fadc[i][j] = new TF1(name,"landau",0,1024);
00134            fadc[i][j]->SetLineColor(1);
00135            fadc[i][j]->SetLineWidth(2);
00136          }
00137      }
00138 
00139 
00140 
00141 
00142 
00143 
00144    /*READ ROOT FILE*/
00145 
00146 
00147    cout<<""<<endl;
00148    cout<<"READ DATA FROM FILE "<<file<<endl;
00149    cout<<""<<endl;
00150 
00151 
00152    TFile *f = new TFile(file);
00153 
00154    TIter nextkey(f.GetListOfKeys());  
00155    TKey *key;
00156         
00157    int nchannel = 0;
00158    int adc = 0;
00159    int hardid = 0;
00160    int chbid = 0;       
00161    int xpos = 0;        
00162    int ypos = 0;        
00163 
00164    int nhot = 0;
00165    int nsuper_hot = 0;
00166    int nhot_tot = 0;
00167    int nsuper_hot_tot = 0;
00168    int ntest[10];
00169    int neff[10];
00170    int nhit[10];
00171    bool hot[10];
00172    bool super_hot[10];
00173    for (int i=0;i<10;i++){
00174      ntest[i] = 0;neff[i] = 0;
00175      nhit[i] = 0;hot[i] = 0;super_hot[i] = 0;}
00176    bool align = 0;
00177    int ntrack = 0;
00178    int strip = 0;
00179 
00180    int xmin = 1000;
00181    int ymin = 1000;
00182 
00183    int xmax = -1;
00184    int ymax = -1;
00185 
00186    while (key = (TKey*)nextkey()) 
00187      {
00188     
00189        TTree *tree = (TTree*)key->ReadObj();                            
00190        MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
00191 
00192        int nEvent = tree.GetEntries();   
00193        cout <<"Number of events = "<<nEvent<<endl;
00194        cout<<endl;
00195 
00196        MTEvent *evt = new MTEvent();
00197        TBranch *branch= tree->GetBranch("MTEvent");
00198        branch->SetAddress(&evt);
00199 
00200        //nEvent = 100;
00201        //for (int i=1;i<nEvent+1;i++)
00202        for (int i=5103;i<5104;i++)
00203          {
00204 
00205            tree.GetEntry(i);
00206            nchannel = evt->GetNchannel();
00207 
00208            if (i%100==0){cout<<"Progress: "<<i/1./nEvent*100<<" \r"<<flush;}
00209 
00210            for (int j=0;j<nchannel;j++)
00211              { 
00212 
00213                MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00214                chbid = channel->GetChamberId();
00215                hardid = channel->GetHardId();
00216 
00217                if (!noisy[chbid-1][hardid])
00218                  {
00219 
00220                    adc = channel->GetAnalogValue();
00221                    hadc[chbid-1][hardid]->Fill(adc);
00222 
00223                    hc_adc[chbid-1]->Fill(hardid,adc);
00224 
00225                    if (chbid<=6 && chbid%2==1)
00226                      {
00227                        strip = channel->GetY()/250 - 0.5 - 272;
00228                        hs_adc[chbid-1]->Fill(strip,adc);
00229                      }
00230                    if (chbid<=6 && chbid%2==0)
00231                      {
00232                        strip = channel->GetX()/250 - 0.5 - 72;
00233                        hs_adc[chbid-1]->Fill(strip,adc);
00234                      }
00235 
00236                    nchan[chbid-1]++;
00237 
00238                    if (adc>adc_cut)
00239                      {
00240                    
00241                        nhit[chbid-1]++;
00242 
00243                        if(rotate90)
00244                          {
00245                            xpos = channel->GetY();
00246                            ypos = 10-channel->GetX();
00247                         }
00248                         else
00249                         {
00250                            xpos = channel->GetX();
00251                            ypos = channel->GetY();
00252                         }
00253 
00254                        if (chbid==7 || chbid==8 || chbid==9)
00255                          {
00256                            hxy_pad_temp->Fill(ypos,xpos);
00257                          }
00258 
00259                        hxy_hit[chbid-1]->Fill(ypos+5,xpos);
00260                        hxy_adc[chbid-1]->Fill(ypos,xpos,adc);
00261                      }
00262                  }
00263              }
00264 
00265 
00266 
00267            /*Find tracks in pad telescope*/
00268 
00269            for (int k=6;k<9;k++)
00270              {
00271 
00272                if (nhit[k]!=0 && nhit[k]<5)
00273                  {
00274                    hot[k] = 1;
00275                    hmult[k]->Fill(nhit[k]);
00276 
00277                    if (nhit[k] == 1)
00278                      {
00279                        super_hot[k] = 1;
00280                      }
00281                  }
00282 
00283                if (hot[k])
00284                  {
00285                    nhot++;
00286                    nhot_tot++;
00287                  }
00288 
00289                if (super_hot[k])
00290                  {
00291                    nsuper_hot++;
00292                    nsuper_hot_tot++;
00293                  }
00294              }
00295 
00296            hnhot->Fill(nhot);
00297            hnsuper_hot->Fill(nsuper_hot);
00298 
00299            align = (hxy_pad_temp->GetRMS(1) == 0 && hxy_pad_temp->GetRMS(2) == 0);
00300            hxy_pad_temp->Reset();
00301 
00302 
00303            if (nsuper_hot>=2 && align)
00304              {
00305                ntrack++;
00306 
00307                if (super_hot[6] && super_hot[7])
00308                  {
00309                    ntest[8]++;
00310                    if (super_hot[8]){neff[8]++;}}
00311 
00312                if (super_hot[6] && super_hot[8])
00313                  {
00314                    ntest[7]++;
00315                    if (super_hot[7]){neff[7]++;}}
00316 
00317                if (super_hot[7] && super_hot[8])
00318                  {
00319                    ntest[6]++;
00320                    if (super_hot[6]){neff[6]++;}}
00321              }
00322 
00323 
00324            nhot = 0;
00325            nsuper_hot = 0;
00326 
00327            for (int k=0;k<9;k++)
00328              {
00329 
00330                hnhit[k]->Fill(nhit[k]);
00331                hnchan[k]->Fill(nchan[k]);
00332                nchan[k] = 0;
00333 
00334                nhit[k] = 0;
00335                hot[k] = 0;
00336                super_hot[k] = 0;
00337 
00338              }
00339          }
00340      }
00341 
00342 
00343    float eff[9];
00344    float eff_err[9];
00345 
00346    for (int k=6;k<9;k++)
00347      {
00348        eff[k] = 100. * neff[k] / ntest[k];
00349        eff_err[k] = 100. * sqrt(eff[k]*0.01 * (1 - eff[k]*0.01) / (1.*ntest[k]));
00350        cout<<eff[k]<<" "<<eff_err[k]<<endl;
00351      }
00352 
00353 
00354 
00355    cout<<endl;
00356    cout<<"N hot events = "<<nhot_tot<<endl;
00357    cout<<"N super hot events = "<<nsuper_hot_tot<<endl;
00358    cout<<"N super hot align events (Ntrack) = "<<ntrack<<endl;
00359    cout<<endl;
00360    cout<<"Efficiency for tracks"<<endl;
00361    cout<<"Chab 7 = "<<neff[6]<<" / "<<ntest[6]<<" = ["<<eff[6]<<" +/- "<<eff_err[6]<<"] %"<<endl;
00362    cout<<"Chab 8 = "<<neff[7]<<" / "<<ntest[7]<<" = ["<<eff[7]<<" +/- "<<eff_err[7]<<"] %"<<endl;
00363    cout<<"Chab 9 = "<<neff[8]<<" / "<<ntest[8]<<" = ["<<eff[8]<<" +/- "<<eff_err[8]<<"] %"<<endl;
00364    cout<<endl;
00365    cout<<"Multiplicity"<<endl;
00366    cout<<"Chab 7 = "<<hmult[6]->GetMean()<<endl;
00367    cout<<"Chab 8 = "<<hmult[7]->GetMean()<<endl;
00368    cout<<"Chab 9 = "<<hmult[8]->GetMean()<<endl;
00369    cout<<endl;
00370 
00371 
00372 
00373  }
00374 
00375 
00376    /*FIT PROFILE DISTRIBUTIONS*/
00377 
00378    if (fit)
00379      {
00380 
00381        cout<<""<<endl;
00382        cout<<"FIT HIT PROFILE"<<endl;
00383        cout<<""<<endl;
00384 
00385        for (int i=0;i<9;i++)
00386          {
00387 
00388            cout<<"******* Chamber "<<i+1<<" *******"<<endl;
00389            cout<<endl;
00390 
00391            if (hxy_hit[i]->GetEntries()!=0)
00392 
00393              {
00394 
00395                fx[i]->SetParameters(hxy_hit[i]->ProjectionX()->GetMaximum(),hxy_hit[i]->ProjectionX()->GetMaximumBin(),hxy_hit[i]->ProjectionX()->GetRMS());
00396                fy[i]->SetParameters(hxy_hit[i]->ProjectionY()->GetMaximum(),hxy_hit[i]->ProjectionY()->GetMaximumBin(),hxy_hit[i]->ProjectionY()->GetRMS());
00397 
00398                hxy_hit[i]->ProjectionX()->Fit(fx[i],"q");
00399                hxy_hit[i]->ProjectionY()->Fit(fy[i],"q");
00400 
00401                cout<<"Horizontal profile (along y)"<<endl;
00402                cout<<"  Mean  = ["<<fx[i]->GetParameter(1)<<" +/- "<<fx[i]->GetParError(1)<<"] cm"<<endl;
00403                cout<<"  Sigma = ["<<fx[i]->GetParameter(2)<<" +/- "<<fx[i]->GetParError(2)<<"] cm"<<endl;
00404 
00405                cout<<"Vertical profile (along x)"<<endl;
00406                cout<<"  Mean  = ["<<fy[i]->GetParameter(1)<<" +/- "<<fy[i]->GetParError(1)<<"] cm"<<endl;
00407                cout<<"  Sigma = ["<<fy[i]->GetParameter(2)<<" +/- "<<fy[i]->GetParError(2)<<"] cm"<<endl;
00408                cout<<endl;
00409 
00410              }
00411 
00412            else
00413              
00414              {
00415 
00416                cout<<"No hits (ADC > "<<adc_cut<<")"<<endl;
00417                cout<<endl;
00418 
00419              }
00420 
00421          }
00422 
00423      }
00424 
00425 
00426 
00427 
00428 
00429 
00430 
00431    if (print)
00432      {
00433 
00434        cout<<endl;
00435        cout<<"PRINT PLOTS"<<endl;
00436        cout<<endl;
00437 
00438 
00439        /*Strip chambers*/
00440 
00441        TCanvas * c1 = new TCanvas("c1","Hit profile telescope - strip chambers",10,10,1500,1000);       
00442        TCanvas * c2 = new TCanvas("c2","Hit profile telescope - strip chambers",10,10,1500,1000);       
00443        TCanvas * c3 = new TCanvas("c3","Hit profile telescope - strip chambers",10,10,1500,1000);
00444        c1->Divide(2,3);
00445        c2->Divide(2,3);
00446        c3->Divide(2,3);
00447 
00448        for (int i=0;i<3;i++)
00449          {
00450            c1->cd(2*i+1);
00451            hxy_hit[2*i]->Draw("zcol");
00452            c1->cd(2*i+2);
00453            hxy_hit[2*i+1]->Draw("zcol");
00454          }
00455 
00456        c1->Print("Plots/hxy_nhit_strip_chb_123456.png");
00457 
00458 
00459        for (int i=0;i<3;i++)
00460          {
00461            c2->cd(2*i+1);
00462            hnchan[2*i]->Draw();
00463            c2->cd(2*i+2);
00464            hnchan[2*i+1]->Draw();
00465          }
00466 
00467        c2->Print("Plots/hnchan_strip_chb_123456.png");
00468 
00469 
00470        for (int i=0;i<3;i++)
00471          {
00472            c3->cd(2*i+1);
00473            hnhit[2*i]->Draw();
00474            c3->cd(2*i+2);
00475            hnhit[2*i+1]->Draw();
00476          }
00477 
00478        c3->Print("Plots/hnhit_strip_chb_123456.png");
00479 
00480 
00481        /*Pad chambers*/
00482 
00483        TCanvas * c4 = new TCanvas("c4","Hit profile telescope - pad chambers",10,10,1000,1000);
00484        TCanvas * c5 = new TCanvas("c5","Hit profile telescope - pad chambers",10,10,1000,1000);
00485        TCanvas * c6 = new TCanvas("c6","Hit profile telescope - pad chambers",10,10,1000,1000);
00486        c4->Divide(1,3);
00487        c5->Divide(1,3);
00488        c6->Divide(1,3);
00489 
00490        for (int i=0;i<3;i++)
00491          {
00492            c4->cd(i+1);
00493            hxy_hit[6+i]->Draw("zcol");
00494          }
00495 
00496        c4->Print("Plots/hxy_nhit_pad_chb_789.png");
00497 
00498 
00499        for (int i=0;i<3;i++)
00500          {
00501            c5->cd(i+1);
00502            hnchan[6+i]->Draw();
00503          }
00504 
00505        c5->Print("Plots/hnchan_pad_chb_789.png");
00506 
00507 
00508        for (int i=0;i<3;i++)
00509          {
00510            c6->cd(i+1);
00511            hnhit[6+i]->Draw();
00512          }
00513 
00514        c6->Print("Plots/hnhit_pad_chb_789.png");
00515 
00516      }
00517 
00518 
00519 
00520    
00521 
00522  }
00523 
00524 
00525
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/bookkeeping/Telescope/Beam/telescope_beam_profile.C
