Micromegas LAPP reconstruction

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