Micromegas LAPP reconstruction

00001     {
00002    gROOT->Reset();
00003 
00004    gStyle->SetPalette(1);
00005    gSystem->Load("/LC/Detecteurs/MicroMegas/Offline/micromegasFrameWork/trunk/lib/libMicro.so");
00006 
00007 
00008    /*m2 prototype time cut*/
00009    int dt_min = 6;
00010    int dt_max = 8;
00011 
00012    /*telescope energy cut*/
00013    int adc_cut = 30;
00014 
00015 
00016 
00017    //READ INFO FROM XML
00018    std::ifstream xml;
00019    xml.open("complete_detector.xml",ifstream::in);
00020    std::string line = "";
00021    std::string lineb = "";
00022    while(line.find("path=")==string::npos || lineb.find("<output")==string::npos)
00023      {
00024         lineb=line;
00025         xml>>line;
00026      }
00027    int start = line.find("\"")+1;
00028    int end =  line.rfind("\"")-6;
00029    TString file = line.substr(start,end);
00030 
00031    for (int i=0;i<5;i++){file.Chop();}
00032    file += "_MGT.root";
00033 
00034 
00035 
00036 
00037    bool noisy[9][32][32];
00038    for (int i=0;i<32;i++){
00039      for (int j=0;j<32;j++){
00040        for (int k=0;k<9;k++){
00041          noisy[k][i][j] = 0;}}}
00042 
00043    /*HV pad*/
00044    noisy[6][5][5] = 1;
00045    noisy[7][5][5] = 1;
00046    noisy[8][5][5] = 1;
00047 
00048    /*floating pad*/
00049    noisy[6][17][5] = 1;
00050    noisy[8][15][9] = 1;
00051 
00052 
00053 
00054    int nchannel = 0;
00055    int hardid = 0;
00056    int chipid = 0;
00057    int nhit = 0;
00058    int nhit_cut = 0;
00059    int neff = 0;
00060    int xpos = 0;  
00061    int ypos = 0;  
00062    int zpos = 0;  
00063    int adc = 0;
00064    int chbid = 0;
00065    int ntree = 0;
00066 
00067    Long64_t t,t0,t1,t2,t3,dt;
00068    
00069    TH2I * hxy_m2_dac0 = new TH2I("hxy_m2_dac0","Hit position distribution;y (pad);x (pad)",96,0,96,96,0,96);
00070    TH2I * hxy_m2_cut_dac0 = new TH2I("hxy_m2_cut_dac0","Hit position distribution - time cut;y (pad);x (pad)",96,0,96,96,0,96);
00071 
00072    TH2I * hxy_m2_dac1 = new TH2I("hxy_m2_dac1","Hit position distribution;y (pad);x (pad)",96,0,96,96,0,96);
00073    TH2I * hxy_m2_cut_dac1 = new TH2I("hxy_m2_cut_dac1","Hit position distribution - time cut;y (pad);x (pad)",96,0,96,96,0,96);
00074 
00075    TH2I * hxy_m2_dac2 = new TH2I("hxy_m2_dac2","Hit position distribution;y (pad);x (pad)",96,0,96,96,0,96);
00076    TH2I * hxy_m2_cut_dac2 = new TH2I("hxy_m2_cut_dac2","Hit position distribution - time cut;y (pad);x (pad)",96,0,96,96,0,96);
00077 
00078    TH1I * hnhit = new TH1I("hnhit","Number of hit distribution;Nhits",10000,0,10000);
00079    TH1I * hnhit_cut = new TH1I("hnhit_cut","Number of hit distribution - time cut;Nhits",10000,0,10000);
00080    TH1I * hdt = new TH1I("hdt","Time to readout disitribution;time (200 ns clock cycle)",200,0,200);
00081    TH2I * hxy_tel[9];
00082    TString title,name;
00083 
00084    for (int i=0;i<9;i++)
00085      {
00086        name = Form("hxy_tel_chamber_%i",i+1);
00087        title = Form("Hit pattern (adc>%i) chamber %i;y (pad),x (pad)",adc_cut,i+1);
00088        hxy_tel[i] = new TH2I(name,title,32,0,32,32,0,32);
00089      }
00090 
00091 
00092 
00093    
00094    TFile *f = new TFile(file);
00095    cout<<endl;
00096    cout<<"Reading file : "<<file<<endl;
00097    cout<<endl;
00098 
00099    TIter nextkey(f.GetListOfKeys());
00100    TKey *key;
00101 
00102 
00103 
00104 
00105    while (key = (TKey*)nextkey())
00106      {
00107 
00108        TTree *tree = (TTree*)key->ReadObj();                
00109        MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
00110 
00111        int nEvent = tree.GetEntries();
00112        ntree++;
00113        cout<<"  Tree "<<ntree<<" contains "<<nEvent<<" events"<<endl;
00114 
00115        MTEvent *evt = new MTEvent();
00116 
00117        TBranch *branch= tree->GetBranch("MTEvent");
00118        branch->SetAddress(&evt);
00119 
00120        for (int i=0;i<nEvent-1;i++)
00121          {
00122 
00123            if (i%100==0){cout<<"Progress: "<<i/1./nEvent*100<<" \r"<<flush;}
00124 
00125            tree.GetEntry(i);
00126            nchannel = evt->GetNchannel();
00127            hnhit->Fill(nchannel);
00128            if (i>0){hnhit_cut->Fill(nhit_cut);nhit_cut = 0;}
00129 
00130            for (int j=0;j<nchannel;j++)
00131              {
00132 
00133                MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00134 
00135                chbid = channel->GetChamberId();
00136                xpos = channel->GetX()*1e-4;
00137                ypos = channel->GetY()*1e-4;
00138 
00139                /*if m2*/
00140                if (chbid == 10)
00141                  {
00142                    t2 = channel->GetBcId_Dif();
00143                    t3 = channel->GetBcId_Hit();
00144                    dt = t2 - t3;
00145                    hdt->Fill(dt);
00146 
00147                    int digit = channel->GetDigitalValue();
00148                    if (digit>=1)
00149                      {
00150                        hxy_m2_dac0->Fill(ypos,xpos);
00151                        if (digit>=2)
00152                          {
00153                            hxy_m2_dac1->Fill(ypos,xpos);
00154                            if (digit==3)
00155                              {
00156                                hxy_m2_dac2->Fill(ypos,xpos);
00157                              }
00158                          }
00159                      }
00160 
00161          
00162                    if (dt>=dt_min && dt<dt_max)
00163                      {
00164                        nhit_cut++;
00165 
00166                        if (digit>=1)
00167                          {
00168                            hxy_m2_cut_dac0->Fill(ypos,xpos);
00169                            if (digit>=2)
00170                              {
00171                                hxy_m2_cut_dac1->Fill(ypos,xpos);
00172                                if (digit==3)
00173                                  {
00174                                    hxy_m2_cut_dac2->Fill(ypos,xpos);
00175                                  }
00176                              }
00177                          }
00178 
00179                        /*estimate of efficiency*/
00180                        if (nhit_cut == 1)
00181                          {
00182 
00183                            neff++;
00184 
00185                          }
00186 
00187                      }
00188                  }
00189                
00190                /*telescope*/
00191                else
00192                  {
00193 
00194                    if (!noisy[chbid-1][ypos][xpos])
00195                      {
00196 
00197                        adc = channel->GetAnalogValue();
00198                        
00199                        if (adc>=adc_cut)
00200                          {
00201                            hxy_tel[chbid-1]->Fill(ypos,xpos);
00202                          }
00203                      }
00204                  }
00205              }
00206          }
00207      }
00208     
00209 
00210 
00211    cout<<endl;
00212    cout<<"Nhot_m2 / Nevent = "<<neff<<" / "<<nEvent<<" = "<<neff*100./nEvent<<" %"<<endl;
00213    cout<<endl;
00214 
00215 
00216 
00217 
00218    int xmax_tel[9];
00219    int ymax_tel[9];
00220    int xmax_m2 = 0;
00221    int ymax_m2 = 0;
00222    int zmax = 0;
00223 
00224 
00225    hxy_m2_cut_dac0->GetMaximumBin(ymax_m2,xmax_m2,zmax);
00226 
00227 
00228    cout<<endl;
00229    cout<<endl;
00230    cout<<"Maximum of hit profile in m2 prototype"<<endl;
00231    cout<<endl;
00232    cout<<"    y-horiz. (pad) = "<<ymax_m2-1<<endl;
00233    cout<<"    x-verti. (pad) = "<<xmax_m2-1<<endl;
00234    cout<<endl;
00235    cout<<endl;
00236 
00237    cout<<"Maximum of hit profile in telescope pad chambers"<<endl;
00238    cout<<endl;
00239 
00240    for (int i=6;i<9;i++)
00241      {
00242  
00243        hxy_tel[i]->GetMaximumBin(ymax_tel[i],xmax_tel[i],zmax);
00244 
00245        cout<<"  Chamber "<<i+1<<endl;
00246        cout<<"    y-horiz. (pad) = "<<ymax_tel[i]-1<<endl;
00247        cout<<"    x-verti. (pad) = "<<xmax_tel[i]-1<<endl;
00248        cout<<endl;
00249 
00250      }
00251 
00252    cout<<endl;
00253 
00254 
00255 
00256     }
00257 
00258 
00259 
00260 
00261 
00262    TCanvas * c0 = new TCanvas("c0","Prototype m2",10,10,1200,800);
00263    c0->Divide(3,2);
00264 
00265    c0->cd(1);
00266    hxy_m2_dac0->Draw("zcol");
00267 
00268    c0->cd(2);
00269    hdt->Draw("zcol");
00270 
00271    c0->cd(3);
00272    hxy_m2_cut_dac0->Draw("zcol");
00273 
00274    c0->cd(4);
00275    hxy_tel[6]->Draw("zcol");
00276 
00277    c0->cd(5);
00278    hxy_tel[7]->Draw("zcol");
00279 
00280    c0->cd(6);
00281    hxy_tel[8]->Draw("zcol");
00282 
00283    c0->Print("Plots/complete_detector_beam_profile_2d.png");
00284 
00285 
00286    TCanvas * c1 = new TCanvas("c1","Prototype m2",10,10,1200,600);
00287    c1->Divide(4,2);
00288 
00289    c1->cd(1);
00290    hxy_m2_cut_dac0->ProjectionX()->SetFillColor(2);
00291    hxy_m2_cut_dac0->ProjectionX()->SetLineColor(1);
00292    hxy_m2_cut_dac0->ProjectionX()->Draw();
00293 
00294    c1->cd(1+4);
00295    hxy_m2_cut_dac0->ProjectionY()->SetFillColor(3);
00296    hxy_m2_cut_dac0->ProjectionY()->SetLineColor(1);
00297    hxy_m2_cut_dac0->ProjectionY()->Draw();
00298 
00299    c1->cd(2);
00300    hxy_tel[6]->ProjectionX()->Draw();
00301    hxy_tel[7]->ProjectionX()->SetLineColor(1);
00302    hxy_tel[6]->ProjectionX()->SetFillColor(2);
00303 
00304    c1->cd(2+4);
00305    hxy_tel[6]->ProjectionY()->Draw();
00306    hxy_tel[6]->ProjectionY()->SetLineColor(1);
00307    hxy_tel[6]->ProjectionY()->SetFillColor(3);
00308 
00309    c1->cd(3);
00310    hxy_tel[7]->ProjectionX()->Draw();
00311    hxy_tel[7]->ProjectionX()->SetLineColor(1);
00312    hxy_tel[7]->ProjectionX()->SetFillColor(2);
00313 
00314    c1->cd(3+4);
00315    hxy_tel[7]->ProjectionY()->Draw();
00316    hxy_tel[7]->ProjectionY()->SetLineColor(1);
00317    hxy_tel[7]->ProjectionY()->SetFillColor(3);
00318 
00319    c1->cd(4);
00320    hxy_tel[8]->ProjectionX()->Draw();
00321    hxy_tel[8]->ProjectionX()->SetLineColor(1);
00322    hxy_tel[8]->ProjectionX()->SetFillColor(2);
00323 
00324    c1->cd(4+4);
00325    hxy_tel[8]->ProjectionY()->Draw();
00326    hxy_tel[8]->ProjectionY()->SetLineColor(1);
00327    hxy_tel[8]->ProjectionY()->SetFillColor(3);
00328 
00329 
00330    c1->Print("Plots/complete_detector_beam_profile_1d.png");
00331 
00332 
00333    cout<<endl;}
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/bookkeeping/Complete_detector/Beam3/complete_detector_profiles.C
