Micromegas LAPP reconstruction

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