Micromegas LAPP reconstruction

00001 {gROOT->Reset();
00002   
00003   gSystem->Load("libMicro.so");
00004   
00005 
00006 
00007   /*Total number of events is 62046*/
00008   /*Number of event with Nchannel != 672 is 17888*/
00009   /*Number of good event (Nch=672) with one single hit (ADC>35) per chamber is 3570*/
00010 
00011 
00012 
00013   bool display = 0;
00014 
00015   
00016   TString path = "/lapp_data/LC/Detecteurs/MicroMegas/data/TB2009/May/Rootfiles/dev/";
00017   //TString path = "/home1/chefdevi/micromegasFrameWork/trunk/src/analyse/Gassiplex/Beam_profile/";
00018   //TString fileName= "beam_profile_alice_pm.root";
00019   TString fileName= "electron_2gev_4absorber.root";
00020   TString fileName2= "histo_beam_profile_alice_pm.root";
00021   
00022   TFile *f = new TFile(path+fileName);
00023   cout<<"reading File : "<<path<<fileName<<endl;
00024   TTree *tree;
00025   tree  = (TTree*) f.Get("t1");cout<<"tree built\n";
00026   int nEvent = tree.GetEntries();
00027   MTEvent *evt = new MTEvent();
00028   TBranch *branch= tree->GetBranch("MTEvent");
00029   branch->SetAddress(&evt);
00030   cout <<"nEvent="<<nEvent<<endl;
00031   int NC = 672;
00032 
00033 
00034 
00035   /*Save histo in a separate ROOT file*/
00036   //TFile *f2 = new TFile(fileName2,"RECREATE");
00037   //f2->Write();
00038   //f2->Close();
00039 
00040 
00041 
00042   /*Distributions of the number of hits*/
00043 
00044   TH1F * hnhit1 = new TH1F("hnhit1","",96,0,96);
00045   TH1F * hnhit2 = new TH1F("hnhit2","",96,0,96);
00046   TH1F * hnhit3 = new TH1F("hnhit3","",96,0,96);
00047   TH1F * hnhit4 = new TH1F("hnhit4","",96*4,0,96*4);
00048 
00049   TH1F * hnhit1_clean = new TH1F("hnhit1_clean","",96*4,0,96*4);
00050   TH1F * hnhit2_clean = new TH1F("hnhit2_clean","",96*4,0,96*4);
00051   TH1F * hnhit3_clean = new TH1F("hnhit3_clean","",96*4,0,96*4);
00052   TH1F * hnhit4_clean = new TH1F("hnhit4_clean","",96*4,0,96*4);
00053 
00054 
00055   /*Distributions of the energy (or ADC)*/
00056 
00057   TH1F * hadc1 = new TH1F("hadc1","",100,0,1024*5);
00058   TH1F * hadc2 = new TH1F("hadc2","",100,0,1024*5);
00059   TH1F * hadc3 = new TH1F("hadc3","",100,0,1024*5);
00060   TH1F * hadc4 = new TH1F("hadc4","",100,0,1024*10);
00061 
00062   TH1F * hadc1_clean = new TH1F("hadc1_clean","",100,0,1024*5);
00063   TH1F * hadc2_clean = new TH1F("hadc2_clean","",100,0,1024*5);
00064   TH1F * hadc3_clean = new TH1F("hadc3_clean","",100,0,1024*5);
00065   TH1F * hadc4_clean = new TH1F("hadc4_clean","",100,0,1024*10);
00066 
00067 
00068   /*Distributions of the hit positions in x and y-planes*/
00069 
00070   TH1F * hxpos1 = new TH1F("hxpos1","",40,0,40);
00071   TH1F * hxpos2 = new TH1F("hxpos2","",40,0,40);
00072   TH1F * hxpos3 = new TH1F("hxpos3","",40,0,40);
00073   TH1F * hxpos4 = new TH1F("hxpos4","",40,0,40);
00074 
00075   TH1F * hxpos1_clean = new TH1F("hxpos1_clean","",40,0,40);
00076   TH1F * hxpos2_clean = new TH1F("hxpos2_clean","",40,0,40);
00077   TH1F * hxpos3_clean = new TH1F("hxpos3_clean","",40,0,40);
00078   TH1F * hxpos4_clean = new TH1F("hxpos4_clean","",40,0,40);
00079 
00080   TH1F * hypos1 = new TH1F("hypos1","",40,0,40);
00081   TH1F * hypos2 = new TH1F("hypos2","",40,0,40);
00082   TH1F * hypos3 = new TH1F("hypos3","",40,0,40);
00083   TH1F * hypos4 = new TH1F("hxpos4","",40,0,40);
00084 
00085   TH1F * hypos1_clean = new TH1F("hypos1_clean","",40,0,40);
00086   TH1F * hypos2_clean = new TH1F("hypos2_clean","",40,0,40);
00087   TH1F * hypos3_clean = new TH1F("hypos3_clean","",40,0,40);
00088   TH1F * hypos4_clean = new TH1F("hypos4_clean","",40,0,40);
00089 
00090 
00091   /*Distributions of the number of hits in the xy-plane*/
00092 
00093   TH2I * hxy1 = new TH2I("hxy1","",32,0,32,32,0,32);
00094   TH2I * hxy2 = new TH2I("hxy2","",32,0,32,32,0,32);
00095   TH2I * hxy3 = new TH2I("hxy3","",32,0,32,32,0,32);
00096   TH2I * hxy4 = new TH2I("hxy4","",32,0,32,32,0,32);
00097 
00098   TH2I * hxy1_all = new TH2I("hxy1_all","",32,0,32,32,0,32);
00099   TH2I * hxy2_all = new TH2I("hxy2_all","",32,0,32,32,0,32);
00100   TH2I * hxy3_all = new TH2I("hxy3_all","",32,0,32,32,0,32);
00101   TH2I * hxy4_all = new TH2I("hxy4_all","",32,0,32,32,0,32);
00102 
00103 
00104 
00105   /*Set Line Width*/
00106   hnhit1->SetLineWidth(2);
00107   hnhit2->SetLineWidth(2);
00108   hnhit3->SetLineWidth(2);
00109   hnhit4->SetLineWidth(2);
00110 
00111   hnhit1_clean->SetLineWidth(2);
00112   hnhit2_clean->SetLineWidth(2);
00113   hnhit3_clean->SetLineWidth(2);
00114   hnhit4_clean->SetLineWidth(2);
00115 
00116   hadc1->SetLineWidth(2);
00117   hadc2->SetLineWidth(2);
00118   hadc3->SetLineWidth(2);
00119   hadc4->SetLineWidth(2);
00120 
00121   hadc1_clean->SetLineWidth(2);
00122   hadc2_clean->SetLineWidth(2);
00123   hadc3_clean->SetLineWidth(2);
00124   hadc4_clean->SetLineWidth(2);
00125 
00126   hxpos1->SetLineWidth(2);
00127   hypos1->SetLineWidth(2);
00128   hxpos2->SetLineWidth(2);
00129   hypos2->SetLineWidth(2);
00130   hxpos3->SetLineWidth(2);
00131   hypos3->SetLineWidth(2);
00132   hxpos4->SetLineWidth(2);
00133   hypos4->SetLineWidth(2);
00134 
00135   hxpos1_clean->SetLineWidth(2);
00136   hypos1_clean->SetLineWidth(2);
00137   hxpos2_clean->SetLineWidth(2);
00138   hypos2_clean->SetLineWidth(2);
00139   hxpos3_clean->SetLineWidth(2);
00140   hypos3_clean->SetLineWidth(2);
00141   hxpos4_clean->SetLineWidth(2);
00142   hypos4_clean->SetLineWidth(2);
00143 
00144 
00145   /*Set Line Color*/
00146   hnhit2->SetLineColor(2);
00147   hnhit3->SetLineColor(4);
00148   hnhit4->SetLineColor(8);
00149 
00150   hnhit2_clean->SetLineColor(2);
00151   hnhit3_clean->SetLineColor(4);
00152   hnhit4_clean->SetLineColor(8);
00153 
00154   hadc2->SetLineColor(2);
00155   hadc3->SetLineColor(4);
00156   hadc4->SetLineColor(8);
00157 
00158   hadc2_clean->SetLineColor(2);
00159   hadc3_clean->SetLineColor(4);
00160   hadc4_clean->SetLineColor(8);
00161 
00162   hxpos2->SetLineColor(2);
00163   hypos2->SetLineColor(2);
00164   hxpos3->SetLineColor(4);
00165   hypos3->SetLineColor(4);
00166   hxpos4->SetLineColor(8);
00167   hypos4->SetLineColor(8);
00168 
00169   hxpos2_clean->SetLineColor(2);
00170   hypos2_clean->SetLineColor(2);
00171   hxpos3_clean->SetLineColor(4);
00172   hypos3_clean->SetLineColor(4);
00173   hxpos4_clean->SetLineColor(8);
00174   hypos4_clean->SetLineColor(8);
00175 
00176 
00177   /*Set Line Style*/
00178   hxpos1_clean->SetLineStyle(2);
00179   hypos1_clean->SetLineStyle(2);
00180   hxpos2_clean->SetLineStyle(2);
00181   hypos2_clean->SetLineStyle(2);
00182 
00183   hxpos3_clean->SetLineStyle(2);
00184   hypos3_clean->SetLineStyle(2);
00185   hxpos4_clean->SetLineStyle(2);
00186   hypos4_clean->SetLineStyle(2);
00187 
00188   hnhit1_clean->SetLineStyle(2);
00189   hnhit2_clean->SetLineStyle(2);
00190   hnhit3_clean->SetLineStyle(2);
00191   hnhit4_clean->SetLineStyle(2);
00192 
00193   hadc1_clean->SetLineStyle(2);
00194   hadc2_clean->SetLineStyle(2);
00195   hadc3_clean->SetLineStyle(2);
00196   hadc4_clean->SetLineStyle(2);
00197 
00198 
00199 
00200 
00201   /*Define variables*/
00202 
00203   int xpos1,ypos1,adc1,sum_adc1;
00204   int xpos2,ypos2,adc2,sum_adc2;
00205   int xpos3,ypos3,adc3,sum_adc3;
00206   int xpos4,ypos4,adc4,sum_adc4;
00207 
00208   vector <int> vxpos1,vypos1,vadc1;
00209   vector <int> vxpos2,vypos2,vadc2;
00210   vector <int> vxpos3,vypos3,vadc3;
00211   vector <int> vxpos4,vypos4,vadc4;
00212 
00213   float avgx[3],avgy[3];
00214   float posz1 = 6;
00215   float posz2 = 13;
00216   float posz3 = 19;
00217   float posz4 = 98;
00218   float avgz[3] = {posz1,posz2,posz3};
00219   float avgx1=0,avgx2=0,avgx3=0,avgx4=0;
00220   float avgy1=0,avgy2=0,avgy3=0,avgy4=0;
00221 
00222   int threshold = 35;
00223   float nhit1=0,nhit2=0,nhit3=0,nhit4=0;
00224   int nclean = 0;
00225   bool clean = 0;
00226 
00227   int event = 0;
00228   int crap_event = 0;
00229 
00230 
00231 
00232 
00233   /*Loop over all events*/
00234 
00235   for( int i=0;i<nEvent;i++){
00236     //for( int i=0;i<10;i++){
00237   //for( int i=event;i<event+1;i++){
00238     
00239     if (display){cout<<"Event # "<<i<<endl;}
00240 
00241     tree.GetEntry(i);
00242     NC = evt->GetNchannel();
00243 
00244     if (NC != 672){crap_event++;}
00245     else{
00246 
00247     clean=0;
00248     nhit1=0,nhit2=0,nhit3=0,nhit4=0;
00249     sum_adc1=0;sum_adc2=0;sum_adc3=0;sum_adc4=0;
00250 
00251     vxpos1.clear();vypos1.clear();vadc1.clear();
00252     vxpos2.clear();vypos2.clear();vadc2.clear();
00253     vxpos3.clear();vypos3.clear();vadc3.clear();
00254     vxpos4.clear();vypos4.clear();vadc4.clear();
00255 
00256     hxy1->Reset();
00257     hxy2->Reset();
00258     hxy3->Reset();
00259     hxy4->Reset();
00260 
00261 
00262     /*Chamber beta 2.1 D*/
00263 
00264     for(int j=0;j<96;j++){
00265 
00266       MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00267       xpos1 = channel.GetX();
00268       ypos1 = channel.GetY();
00269       adc1  = channel.GetValue();
00270       
00271       if (adc1>=threshold){
00272         
00273         vxpos1.push_back(xpos1);
00274         vypos1.push_back(ypos1);
00275         vadc1.push_back(adc1);
00276         sum_adc1+=adc1;
00277 
00278         if (display){cout<<" x1="<<xpos1<<"  y1="<<ypos1<<"  adc1="<<adc1<<endl;}
00279 
00280         nhit1++;
00281         hxpos1->Fill(xpos1);
00282         hypos1->Fill(ypos1);
00283         hxy1->SetBinContent(xpos1+1,ypos1+1,hxy1->GetBinContent(xpos1+1,ypos1+1)+1);
00284         hxy1_all->SetBinContent(xpos1+1,ypos1+1,hxy1_all->GetBinContent(xpos1+1,ypos1+1)+1);}}
00285 
00286 
00287     /*Chamber beta 2.1 E*/
00288 
00289     for(int j=0+96;j<96+96;j++){
00290 
00291       MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00292       xpos2 = channel.GetX();
00293       ypos2 = channel.GetY();
00294       adc2  = channel.GetValue();
00295 
00296       if (adc2>=threshold){
00297 
00298         vxpos2.push_back(xpos2);
00299         vypos2.push_back(ypos2);
00300         vadc2.push_back(adc2);
00301         sum_adc2+=adc2;
00302 
00303         if (display){cout<<" x2="<<xpos2<<"  y2="<<ypos2<<"  adc2="<<adc2<<endl;}
00304 
00305         nhit2++;
00306         hxpos2->Fill(xpos2);
00307         hypos2->Fill(ypos2);
00308         hxy2->SetBinContent(xpos2+1,ypos2+1,hxy2->GetBinContent(xpos2+1,ypos2+1)+1);
00309         hxy2_all->SetBinContent(xpos2+1,ypos2+1,hxy2_all->GetBinContent(xpos2+1,ypos2+1)+1);}}
00310     
00311 
00312 
00313     /*Chamber beta 2.1 A*/
00314 
00315     for(int j=0+96+96;j<96+96+96;j++){
00316 
00317       MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00318       xpos3 = channel.GetX();
00319       ypos3 = channel.GetY();
00320       adc3  = channel.GetValue();
00321 
00322       if (adc3>=threshold){
00323 
00324         vxpos3.push_back(xpos3);
00325         vypos3.push_back(ypos3);
00326         vadc3.push_back(adc3);
00327         sum_adc3+=adc3;
00328 
00329         if (display){cout<<" x3="<<xpos3<<"  y3="<<ypos3<<"  adc3="<<adc3<<endl;}
00330 
00331         nhit3++;
00332         hxpos3->Fill(xpos3);
00333         hypos3->Fill(ypos3);
00334         hxy3->SetBinContent(xpos3+1,ypos3+1,hxy3->GetBinContent(xpos3+1,ypos3+1)+1);
00335         hxy3_all->SetBinContent(xpos3+1,ypos3+1,hxy3_all->GetBinContent(xpos3+1,ypos3+1)+1);}}
00336 
00337 
00338     
00339 
00340 
00341     /*chamber beta 2.4*/
00342 
00343     for(int j=287;j<672;j++){
00344       //for(int j=287;j<668;j++){
00345       MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00346       xpos4 = channel.GetX();
00347       ypos4 = channel.GetY();
00348       adc4  = channel.GetValue();
00349 
00350       if (adc4>=threshold){       
00351 
00352         vxpos4.push_back(xpos4);
00353         vypos4.push_back(ypos4);
00354         vadc4.push_back(adc4);
00355         sum_adc4+=adc4;
00356 
00357         if (display){cout<<" x4="<<xpos4<<"  y4="<<ypos4<<"  adc4="<<adc4<<endl;}
00358 
00359         nhit4++;
00360         hxpos4->Fill(xpos4);
00361         hypos4->Fill(ypos4);
00362         hxy4->SetBinContent(xpos4+1,ypos4+1,hxy4->GetBinContent(xpos4+1,ypos4+1)+1);
00363         hxy4_all->SetBinContent(xpos4+1,ypos4+1,hxy4_all->GetBinContent(xpos4+1,ypos4+1)+1);}}
00364 
00365 
00366     
00367     hnhit1->Fill(nhit1);
00368     hnhit2->Fill(nhit2);
00369     hnhit3->Fill(nhit3);
00370     hnhit4->Fill(nhit4);
00371 
00372     hadc1->Fill(sum_adc1);
00373     hadc2->Fill(sum_adc2);
00374     hadc3->Fill(sum_adc3);
00375     hadc4->Fill(sum_adc4);
00376 
00377 
00378 
00379 
00380     /*clean event if one and only one hit in each beta 2.1 chamber*/
00381     clean = (nhit1==1) && (nhit2==1) && (nhit3==1) && (nhit4==1);
00382     
00383     /*clean event if at least one hit in each beta 2.1 chamber*/
00384     //clean = (nhit1>=1) && (nhit2>=1) && (nhit3>=1);
00385 
00386     /*clean event if at least one hit in each beta 2.1 chamber but less than 10*/
00387     //clean = (nhit1>=1) && (nhit2>=1) && (nhit3>=1) && (nhit1<10) && (nhit2<10) && (nhit3<10);        
00388   
00389 
00390     
00391     
00392     if (clean){
00393 
00394       if (display){cout<<"clean event!"<<endl;}
00395 
00396       nclean++;
00397       
00398       avgx1=0;avgx2=0;avgx3=0;avgx4=0;
00399       avgy1=0;avgy2=0;avgy3=0;avgy4=0;
00400 
00401       hnhit1_clean->Fill(nhit1);
00402       hnhit2_clean->Fill(nhit2);
00403       hnhit3_clean->Fill(nhit3);
00404       hnhit4_clean->Fill(nhit4);
00405 
00406       hadc1_clean->Fill(sum_adc1);
00407       hadc2_clean->Fill(sum_adc2);
00408       hadc3_clean->Fill(sum_adc3);
00409       hadc4_clean->Fill(sum_adc4);
00410 
00411       for (int k=0;k<nhit1;k++){
00412         hxpos1_clean->Fill(vxpos1[k]);
00413         hypos1_clean->Fill(vypos1[k]);}
00414 
00415       for (int k=0;k<nhit2;k++){
00416         hxpos2_clean->Fill(vxpos2[k]);
00417         hypos2_clean->Fill(vypos2[k]);}
00418 
00419       for (int k=0;k<nhit3;k++){
00420         hxpos3_clean->Fill(vxpos3[k]);
00421         hypos3_clean->Fill(vypos3[k]);}
00422 
00423       for (int k=0;k<nhit4;k++){
00424         hxpos4_clean->Fill(vxpos4[k]);
00425         hypos4_clean->Fill(vypos4[k]);}}
00426 
00427 
00428     if (display){cout<<nhit1<<" "<<nhit2<<" "<<nhit3<<" "<<nhit4<<endl;}}}
00429 
00430 
00431       
00432 
00433 
00434 
00435   cout<<"ntotal = "<<nEvent<<"  nclean = "<<nclean<<endl;
00436 
00437   hxy4_all->Draw("zcol");
00438 
00439 
00440 
00441 
00442 
00443 
00444 
00445 
00446 
00447   /*Plot distributions*/
00448 
00449 
00450   bool nhit_distribution = 0;
00451   bool adc_distribution = 0;
00452   bool chamber_xy_display = 0;
00453   bool xyz_plot = 0;
00454 
00455   bool number_hit_distrib_4_chamber = 0;
00456   bool number_hit_distrib_beta24 = 0;
00457 
00458   bool x_y_hit_distrib_beta21 = 0;
00459   bool x_y_hit_distrib_beta22 = 0;
00460   bool x_y_hit_distrib_beta23 = 0;
00461   bool x_y_hit_distrib_beta24 = 0;
00462 
00463 
00464   if (nhit_distribution){
00465 
00466     TCanvas * myc = new TCanvas("myc");
00467     myc->Divide(2,2);
00468     myc->cd(1);
00469     hnhit1->Draw();
00470     hnhit1_clean->Draw("same");
00471     myc->cd(2);
00472     hnhit2->Draw();
00473     hnhit2_clean->Draw("same");
00474     myc->cd(3);
00475     hnhit3->Draw();
00476     hnhit3_clean->Draw("same");
00477     myc->cd(4);
00478     hnhit4->Draw();
00479     hnhit4_clean->Draw("same");}
00480 
00481   if (adc_distribution){
00482 
00483     TCanvas * myc = new TCanvas("myc");
00484     myc->Divide(2,2);
00485     myc->cd(1);
00486     hadc1->Draw();
00487     hadc1_clean->Draw("same");
00488     myc->cd(2);
00489     hadc2->Draw();
00490     hadc2_clean->Draw("same");
00491     myc->cd(3);
00492     hadc3->Draw();
00493     hadc3_clean->Draw("same");
00494     myc->cd(4);
00495     hadc4->Draw();
00496     hadc4_clean->Draw("same");}
00497 
00498 
00499   if (chamber_xy_display){
00500 
00501     TCanvas * myc3 = new TCanvas("myc3");
00502     myc3->Divide(2,2);
00503     myc3->cd(1);
00504     hxy1->Draw("zcol");
00505     myc3->cd(2);
00506     hxy2->Draw("zcol");
00507     myc3->cd(3);
00508     hxy3->Draw("zcol");
00509     myc3->cd(4);
00510     hxy4->Draw("zcol");}
00511 
00512 
00513   if (xyz_plot){
00514     
00515     TCanvas * myc4 = new TCanvas("myc4");
00516 
00517     tgxy->SetMarkerStyle(21);
00518     tgzx->SetMarkerStyle(21);
00519     tgzy->SetMarkerStyle(21);
00520     
00521     myc4->Divide(2,2);
00522     myc4->cd(1);
00523     tgxy->Draw("ap");
00524     myc4->cd(2);
00525     tgzx->Draw("ap");
00526     myc4->cd(3);
00527     tgzy->Draw("ap");}
00528 
00529 
00530   if (number_hit_distrib_4_chamber){
00531     
00532     hnhit1->DrawNormalized();
00533     hnhit2->DrawNormalized("same");
00534     hnhit3->DrawNormalized("same");
00535     hnhit4->DrawNormalized("same");}
00536 
00537   if (number_hit_distrib_beta24){
00538 
00539     hnhit4->DrawNormalized();
00540     hnhit4_clean->DrawNormalized("same");}
00541 
00542   if (x_y_hit_distrib_beta21){
00543 
00544     hxpos1->DrawNormalized();
00545     hxpos1_clean->DrawNormalized("same");
00546     hypos1->DrawNormalized("same");
00547     hypos1_clean->DrawNormalized("same");}
00548 
00549   if (x_y_hit_distrib_beta22){
00550 
00551     hxpos2->DrawNormalized();
00552     hxpos2_clean->DrawNormalized("same");
00553     hypos2->DrawNormalized("same");
00554     hypos2_clean->DrawNormalized("same");}
00555 
00556   if (x_y_hit_distrib_beta23){
00557 
00558     hxpos3->DrawNormalized();
00559     hxpos3_clean->DrawNormalized("same");
00560     hypos3->DrawNormalized("same");
00561     hypos3_clean->DrawNormalized("same");}
00562 
00563   if (x_y_hit_distrib_beta24){
00564 
00565     hxpos4_clean->DrawNormalized();
00566     hxpos4->DrawNormalized("same");
00567     hypos4_clean->DrawNormalized("same");
00568     hypos4->DrawNormalized("same");}
00569 
00570 
00571 }
/Users/jacquem/NetBeansProjects/MicromegasFramework/trunk/src/analyse/Beam_profile.C
