Micromegas LAPP reconstruction

00001 {gROOT->Reset();
00002   
00003   gSystem->Load("/home1/chefdevi/micromegasFrameWork/trunk/libMicro.so");
00004   gStyle->SetPalette(1);
00005   //gStyle->SetOptStat(0);
00006 
00007 
00008   bool display = 0;
00009   bool plot = 1;
00010 
00011 
00012 
00013   TString file = "acq_31102009_2128_2253_1";
00014 
00015   TString filename = "Root_files/"+ file +".root";
00016 
00017 
00018 
00019   int xpos,ypos,zpos,val;
00020   Long64_t t,t1,t2,t3,dt;
00021   int nchannel = 0;
00022   bool time_cut = 0;
00023   int hardid = 0;
00024   int softid = 0;
00025   int nhit = 0;
00026 
00027  
00028 
00029   TH1I * hx_1 = new TH1I("hx_1","X distribution of hits in chamber 1",8,0,8);
00030   TH1I * hx_2 = new TH1I("hx_2","X distribution of hits in chamber 2",8,0,8);
00031   TH1I * hx_3 = new TH1I("hx_3","X distribution of hits in chamber 3",8,0,8);
00032   TH1I * hx_4 = new TH1I("hx_4","X distribution of hits in chamber 4",8,0,8);
00033 
00034   TH1I * hy_1 = new TH1I("hy_1","Y distribution of hits in chamber 1",8,0,8);
00035   TH1I * hy_2 = new TH1I("hy_2","Y distribution of hits in chamber 2",8,0,8);
00036   TH1I * hy_3 = new TH1I("hy_3","Y distribution of hits in chamber 3",8,0,8);
00037   TH1I * hy_4 = new TH1I("hy_4","Y distribution of hits in chamber 4",8,0,8);
00038 
00039   TH2I * hxy_1 = new TH2I("hxy_1","XY distribution of hits in chamber 1",8,0,8,8,0,8);
00040   TH2I * hxy_2 = new TH2I("hxy_2","XY distribution of hits in chamber 2",8,0,8,8,0,8);
00041   TH2I * hxy_3 = new TH2I("hxy_3","XY distribution of hits in chamber 3",8,0,8,8,0,8);
00042   TH2I * hxy_4 = new TH2I("hxy_4","XY distribution of hits in chamber 4",8,0,8,8,0,8);
00043 
00044   TH1I * hz = new TH1I("hz","Z distribution of hits",25,-5,20);
00045 
00046   TH1I * hc_1 = new TH1I("hc_1","Channel number of hits in chamber 1",64,0,64);
00047   TH1I * hc_2 = new TH1I("hc_2","Channel number of hits in chamber 2",64,0,64);
00048   TH1I * hc_3 = new TH1I("hc_3","Channel number of hits in chamber 3",64,0,64);
00049   TH1I * hc_4 = new TH1I("hc_4","Channel number of hits in chamber 4",64,0,64);
00050 
00051   int bcid_hit_max = 10000;
00052   TH1I * hbcid_hit_1 = new TH1I("hbcid_hit_1","BCID_hit distribution in chamber 1",1000,0,bcid_hit_max);
00053   TH1I * hbcid_hit_2 = new TH1I("hbcid_hit_2","BCID_hit distribution in chamber 2",1000,0,bcid_hit_max);
00054   TH1I * hbcid_hit_3 = new TH1I("hbcid_hit_3","BCID_hit distribution in chamber 3",1000,0,bcid_hit_max);
00055   TH1I * hbcid_hit_4 = new TH1I("hbcid_hit_4","BCID_hit distribution in chamber 4",1000,0,bcid_hit_max);
00056 
00057   int bcid_dif_max = 2000;
00058   TH1I * hbcid_dif_1 = new TH1I("hbcid_dif_1","BCID_dif distribution in chamber 1",1000,0,bcid_dif_max);
00059   TH1I * hbcid_dif_2 = new TH1I("hbcid_dif_2","BCID_dif distribution in chamber 2",1000,0,bcid_dif_max);
00060   TH1I * hbcid_dif_3 = new TH1I("hbcid_dif_3","BCID_dif distribution in chamber 3",1000,0,bcid_dif_max);
00061   TH1I * hbcid_dif_4 = new TH1I("hbcid_dif_4","BCID_dif distribution in chamber 4",1000,0,bcid_dif_max);
00062 
00063   int bcid_abs_max = 1e12;
00064   TH1I * hbcid_abs_1 = new TH1I("hbcid_abs_1","BCID_abs distribution in chamber 1",1000,0,bcid_abs_max);
00065   TH1I * hbcid_abs_2 = new TH1I("hbcid_abs_2","BCID_abs distribution in chamber 2",1000,0,bcid_abs_max);
00066   TH1I * hbcid_abs_3 = new TH1I("hbcid_abs_3","BCID_abs distribution in chamber 3",1000,0,bcid_abs_max);
00067   TH1I * hbcid_abs_4 = new TH1I("hbcid_abs_4","BCID_abs distribution in chamber 4",1000,0,bcid_abs_max);
00068 
00069   TH1I * hnhit_1 = new TH1I("hnhit_1","Number of hits per event (per DIF readout) in chamber 1",100,0,100);
00070   TH1I * hnhit_2 = new TH1I("hnhit_2","Number of hits per event (per DIF readout) in chamber 2",100,0,100);
00071   TH1I * hnhit_3 = new TH1I("hnhit_3","Number of hits per event (per DIF readout) in chamber 3",100,0,100);
00072   TH1I * hnhit_4 = new TH1I("hnhit_4","Number of hits per event (per DIF readout) in chamber 4",100,0,100);
00073 
00074 
00075 
00076   int n1=0;
00077   int n2=0;
00078   int n3=0;
00079   int n4=0;
00080   int one_hit_each_chamber = 0;
00081 
00082 
00083 
00084   TFile *f = new TFile(filename);
00085   cout<<"reading File : "<<filename<<endl;
00086 
00087 
00088   TIter nextkey(f.GetListOfKeys());
00089   TKey *key;
00090 
00091   while (key = (TKey*)nextkey()) {
00092     
00093     TTree *tree = (TTree*)key->ReadObj();                
00094     MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
00095     run->Info();
00096 
00097 
00098   int nEvent = tree.GetEntries();
00099   cout <<"nEvent="<<nEvent<<endl;
00100   float nevent = nEvent;
00101 
00102 
00103   MTEvent *evt = new MTEvent();
00104 
00105   TBranch *branch= tree->GetBranch("MTEvent");
00106   branch->SetAddress(&evt);
00107 
00108 
00109 
00110   bool cut = true;
00111 
00112   for (int i=0;i<nevent;i++){  
00113     //for (int i=nevent-10;i<nevent;i++){  
00114 
00115     if (i%10 == 0){cout<<i/nevent*100<<" %             "<<"\r"<<flush;}
00116 
00117     tree.GetEntry(i);
00118     nchannel = evt->GetNchannel();
00119     nhit+=nchannel;
00120 
00121     //cout<<"event number = "<<i<<"  Number of hits = "<<nchannel<<endl;
00122 
00123 
00124     n1=0;n2=0;n3=0;n4=0;
00125 
00126 
00127     for (int j=0;j<nchannel;j++){
00128 
00129     MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00130     xpos = channel->GetX();
00131     ypos = channel->GetY();
00132     zpos = channel->GetZ();
00133     val = channel->GetDigitalValue();
00134     hardid = channel->GetHardId();
00135     softid = channel->GetSoftId();
00136 
00137     t1 = channel->GetBcId_Abs();
00138     t2 = channel->GetBcId_Dif();
00139     t3 = channel->GetBcId_Hit();
00140     dt = (t2-t3);
00141     t = (t1 - (t2-t3));
00142 
00143         cut = ((t2-t3)<5 && (t2-t3)>-5 && t3>0 && t3<1630);
00144     if(cut){
00145 
00146     hz->Fill(zpos);
00147 
00148     if (zpos == 0) {
00149       hx_1->Fill(xpos);hy_1->Fill(ypos);hxy_1->Fill(xpos,ypos);hc_1->Fill(hardid);hbcid_hit_1->Fill(t3);hbcid_dif_1->Fill(t2);hbcid_abs_1->Fill(t1);hnhit_1->Fill(nchannel);if (hardid!=11 && hardid!=57){n1++;}}
00150 
00151     if (zpos == 5) {
00152       hx_2->Fill(xpos);hy_2->Fill(ypos);hxy_2->Fill(xpos,ypos);hc_2->Fill(hardid);hbcid_hit_2->Fill(t3);hbcid_dif_2->Fill(t2);hbcid_abs_2->Fill(t1);hnhit_2->Fill(nchannel);if (hardid!=11 && hardid!=57){n2++;}}
00153 
00154     if (zpos == 10){
00155       hx_3->Fill(xpos);hy_3->Fill(ypos);hxy_3->Fill(xpos,ypos);hc_3->Fill(hardid);hbcid_hit_3->Fill(t3);hbcid_dif_3->Fill(t2);hbcid_abs_3->Fill(t1);hnhit_3->Fill(nchannel);if (hardid!=11 && hardid!=57){n3++;}}
00156 
00157     if (zpos == 15){
00158       hx_4->Fill(xpos);hy_4->Fill(ypos);hxy_4->Fill(xpos,ypos);hc_4->Fill(hardid);hbcid_hit_4->Fill(t3);hbcid_dif_4->Fill(t2);hbcid_abs_4->Fill(t1);hnhit_4->Fill(nchannel);if (hardid!=11 && hardid!=57){n4++;}}
00159 
00160     time_cut = (dt>=2 && dt<=8);
00161 
00162     if (display){cout<<" x = "<<xpos<<" y = "<<ypos<<" z = "<<zpos<<" val = "<<val<<" hardid = "<<hardid<<" softid = "<<softid<<" bcid_abs = "<<t1<<" bcid_dif = "<<t2<<" bcid_hit = "<<t3<<endl;}}}}}
00163 
00164 
00165   cout<<"Number of hits = "<<nhit<<endl;
00166 
00167 
00168 
00169 
00170   if (plot){
00171 
00172     TCanvas * c1 = new TCanvas("c1","",10,10,1600,400);
00173     c1->Divide(4,1);
00174     c1->cd(1);
00175     hxy_1->Draw("zcol");
00176     c1->cd(2);
00177     hxy_2->Draw("zcol");
00178     c1->cd(3);
00179     hxy_3->Draw("zcol");
00180     c1->cd(4);
00181     hxy_4->Draw("zcol");
00182 
00183     TCanvas * c2 = new TCanvas("c2","",10,10,1600,400);
00184     c2->Divide(4,1);
00185     c2->cd(1);
00186     hc_1->Draw();
00187     c2->cd(2);
00188     hc_2->Draw();
00189     c2->cd(3);
00190     hc_3->Draw();
00191     c2->cd(4);
00192     hc_4->Draw();
00193 
00194     TCanvas * c3 = new TCanvas("c3","",10,10,1600,400);
00195     c3->Divide(4,1);
00196     c3->cd(1);
00197     hbcid_hit_1->Draw();
00198     c3->cd(2);
00199     hbcid_hit_2->Draw();
00200     c3->cd(3);
00201     hbcid_hit_3->Draw();
00202     c3->cd(4);
00203     hbcid_hit_4->Draw();
00204 
00205     TCanvas * c4 = new TCanvas("c4","",10,10,600,600);
00206     hz->Draw();
00207 
00208     TCanvas * c5 = new TCanvas("c5","",10,10,1600,400);
00209     c5->Divide(4,1);
00210     c5->cd(1);
00211     hnhit_1->Draw();
00212     c5->cd(2);
00213     hnhit_2->Draw();
00214     c5->cd(3);
00215     hnhit_3->Draw();
00216     c5->cd(4);
00217     hnhit_4->Draw();
00218 
00219     TCanvas * c6 = new TCanvas("c6","",10,10,1600,400);
00220     c6->Divide(4,1);
00221     c6->cd(1);
00222     hbcid_dif_1->Draw();
00223     c6->cd(2);
00224     hbcid_dif_2->Draw();
00225     c6->cd(3);
00226     hbcid_dif_3->Draw();
00227     c6->cd(4);
00228     hbcid_dif_4->Draw();
00229 
00230     TCanvas * c7 = new TCanvas("c7","",10,10,1600,400);
00231     c7->Divide(4,1);
00232     c7->cd(1);
00233     hbcid_abs_1->Draw();
00234     c7->cd(2);
00235     hbcid_abs_2->Draw();
00236     c7->cd(3);
00237     hbcid_abs_3->Draw();
00238     c7->cd(4);
00239     hbcid_abs_4->Draw();
00240 
00241 
00242     TString plotname1 = "Plots/"+ file +"_xy_hit.png";
00243     TString plotname2 = "Plots/"+ file +"_channel_number_hit.png";
00244     TString plotname3 = "Plots/"+ file +"_bcid_hit.png";
00245     TString plotname4 = "Plots/"+ file +"_z_hit.png";
00246     TString plotname5 = "Plots/"+ file +"_number_of_hit.png";
00247     TString plotname6 = "Plots/"+ file +"_bcid_dif.png";
00248     TString plotname7 = "Plots/"+ file +"_bcid_abs.png";
00249 
00250     c1->Print(plotname1.Data());
00251     c2->Print(plotname2.Data());
00252     c3->Print(plotname3.Data());
00253     c4->Print(plotname4.Data());
00254     c5->Print(plotname5.Data());
00255     c6->Print(plotname6.Data());
00256     c7->Print(plotname7.Data());
00257 
00258 
00259 
00260   }
00261 
00262 
00263 
00264 
00265 }
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Ambroise/Analysis_dirac_AE.C
