Micromegas LAPP reconstruction

int main	(	int	argc,
		char **	argv
	)
Definition at line 43 of file geomVerif.cpp.
00045 {
00046 
00047   /*m2 prototype time cut*/
00048   int dt_min = 6;
00049   int dt_max = 9;
00050   
00051   /*telescope energy cut*/
00052   int adc_cut = 30;
00053 
00054   if ( argc != 2  && argc != 3) {
00055     FILE_LOG(logERROR)  << "usage:" << endl;
00056     FILE_LOG(logERROR)  << "geomVerif rootFile [stripZRotation(90)]" << endl;
00057     exit(1);
00058   }
00059 
00060 
00061   string file = argv[1];
00062 
00063   bool rotation =  false;
00064   if ( argc == 3)
00065   {
00066     if ( atoi(argv[2]) == 90 )
00067     {rotation = true; }
00068     else { cout << " stripZRotation should be 90 or NULL "  << endl; exit(-1);}
00069   }
00070     
00071 
00072   TApplication *theApp;
00073 
00074   theApp = new TRint("App", &argc, argv);
00075 
00076   bool noisy[9][96];
00077   for (int i=0;i<9;i++){ for (int j=0;j<96;j++){ noisy[i][j] = 0;}}
00078 
00079   /*HV pad*/
00080   noisy[6][77] = 1;
00081   noisy[7][77] = 1;
00082   noisy[8][77] = 1;
00083 
00084   /*floating pad*/
00085   noisy[6][9] = 1;
00086   noisy[8][33] = 1;
00087 
00088 
00089   bool super_hot[10];
00090   int nsuper_hot_tel = 0;
00091   int nhit[10] = {0};
00092   bool align = 0;
00093   bool clusterInStrip[2] = {false,false};
00094 
00095   int nchannel = 0;
00096   int hardid = 0;
00097   int chipid = 0;
00098 
00099   Float_t xpos = 0.;  
00100   Float_t ypos = 0.;  
00101   Float_t zpos = 0.;  
00102 
00103   int adc = 0;
00104   int chbid = 0;
00105   int ntree = 0;
00106 
00107   Float_t pad_x = 0.;
00108   Float_t pad_y = 0.;
00109 
00110   Float_t xStripNum = 0.;  
00111   Float_t yStripNum = 0.;  
00112 
00113   Int_t pad[2][MAXSTRIPENTRY];
00114   Int_t strip[2][MAXSTRIPENTRY];
00115   Int_t nEntries[2] = {0,0};
00116 
00117   Long64_t t,t0,t1,t2,t3,dt;
00118 
00119   TH2I * hxy_tel_pad_temp = new TH2I("hxy_tel_pad_temp","Hit position distribution in telescope pad chamber for super hot events (Nhit=1);y (pad);x (pad)",32,0,32,32,0,32);
00120   TH2I * hx_cor = new TH2I("hx_cor","Correlation Prototype - telescope;Xtarget (pad);Xm2 (pad)",100,0,1000000,100,0,1000000);
00121   TH2I * hy_cor = new TH2I("hy_cor","Correlation Prototype - telescope;Ytarget (pad);Ym2 (pad)",100,0,1000000,100,0,1000000);
00122   TH1I * hx_res = new TH1I("hx_res","Vertical residuals Prototype - telescope;Xm2 - Xtarget (pad)",101,-500000,500000);
00123   TH1I * hy_res = new TH1I("hy_res","Horizontal residuals Prototype - telescope;Ym2 - Ytarget (pad)",101,-500000,500000);
00124   
00125  
00126   TFile *f = new TFile(file.c_str());
00127   cout<<endl;
00128   cout<<"Reading file : "<<file<<endl;
00129   cout<<endl;
00130 
00131   if (f == NULL) { exit(-1); }
00132   TIter nextkey(f->GetListOfKeys());
00133   TKey *key;
00134 
00135   while (key = (TKey*)nextkey())
00136   {
00137 
00138     TTree *tree = (TTree*)key->ReadObj();                
00139     if ( tree == NULL ) { exit(-1); }
00140     MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
00141 
00142     int nEvent = tree->GetEntries();
00143     ntree++;
00144     cout<<"  Tree "<<ntree<<" contains "<<nEvent<<" events"<<endl;
00145 
00146     MTEvent *evt = new MTEvent();
00147 
00148     TBranch *branch= tree->GetBranch("MTEvent");
00149     branch->SetAddress(&evt);
00150 
00151     for (int i=0;i<nEvent-1 ;i++)
00152     {
00153       if (i%100==0){cout<<"Progress: "<<i/1./nEvent*100<<" \r"<<flush;}
00154       set<Float_t> stripX;
00155       set<Float_t> stripY;
00156       clusterInStrip[X] = false;
00157       clusterInStrip[Y] = false;
00158 
00159       tree->GetEntry(i);
00160       nchannel = evt->GetNchannel();
00161 
00162       nsuper_hot_tel = 0;
00163 
00164       for (int k=0;k<10;k++){nhit[k] = 0; super_hot[k] = 0;}
00165 
00166       // Pour chaque evenement:
00167       // 1 On cherche dans le telescope 1 hit par chambre
00168       // Avec iun filtre sutr les chambres brillante
00169       // et une coupuresur la valeur de l'adc
00170 
00171       /*Look in telescope*/
00172       for (int j=0;j<nchannel;j++)
00173       {
00174         MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00175         chbid = channel->GetChamberId();
00176         if (chbid < 10)
00177         {
00178           if (!noisy[chbid-1][hardid])
00179           {
00180             xpos = channel->GetX();//*1e-4;
00181             ypos = channel->GetY();//*1e-4;
00182             adc = channel->GetAnalogValue();
00183             if (adc>=adc_cut)
00184             {              
00185               nhit[chbid-1]++;
00186               if (chbid == 7 || chbid == 8 || chbid == 9)
00187               {
00188                 hxy_tel_pad_temp->Fill(ypos,xpos);
00189               }
00190               if (chbid == 8)
00191               {
00192                 pad_x = xpos;
00193                 pad_y = ypos;
00194               }
00195               if (  (  chbid == 1 && rotation == false) || ( chbid == 2 && rotation == true)) // Strip gassiplex HORIZONTAL
00196 
00197               {
00198                  stripX.insert(xpos); 
00199               }
00200               if (  (  chbid == 2 && rotation == false) || ( chbid == 1 && rotation == true)) // Strip gassiplex HORIZONTAL
00201               {
00202                  stripY.insert(ypos); 
00203               }
00204             } // end if (adc>=adc_cut)
00205           }  // end if (!noisy[chbid-1][hardid])
00206         } // end  if (chbid < 10)
00207       } // end for (int j=0;j<nchannel;j++)
00208 
00209       // Find cluster in strip chamber
00210       // eliminate noise by filter on particle cluster formed with 3 to 5 consecutive strips
00211       if ( stripX.size() >= 3 &&  stripX.size() <=5 )
00212       {
00213   
00214 // Over X
00215         clusterInStrip[X] = true;
00216         Float_t lastStrip = -1.;
00217         xStripNum = 0.;
00218 
00219         for( set<Float_t>::iterator ii=stripX.begin(); ii!=stripX.end(); ++ii)
00220         {
00221           Float_t foo = (*ii ) / STRIPWIDTH;
00222           xStripNum = xStripNum + foo;
00223           if ( lastStrip != -1 && foo != lastStrip + 1. )
00224           {
00225             clusterInStrip[X] = false; 
00226           }
00227           lastStrip = foo;
00228           
00229         }
00230         xStripNum = xStripNum / stripX.size();
00231       }
00232 
00233       if ( stripY.size() >= 3 &&  stripY.size() <=5 )
00234       {
00235 
00236 // Over Y
00237         clusterInStrip[Y] = true;
00238         Float_t lastStrip = -1.;
00239         yStripNum = 0.;
00240 
00241         for( set<Float_t>::iterator ii=stripY.begin(); ii!=stripY.end(); ++ii)
00242         {
00243           Float_t foo = (*ii ) / STRIPWIDTH;
00244           yStripNum = yStripNum + foo;
00245           if ( lastStrip != -1 && foo != lastStrip + 1. )
00246           {
00247             clusterInStrip[Y] = false;
00248           }
00249           lastStrip = foo;
00250 
00251         }
00252         yStripNum = yStripNum / stripY.size();
00253  
00254 
00255       }
00256  
00257 
00258       //Find tracks in pad telescope
00259       //Only tracks which have only one Pad hit per chamber  
00260       for (int k=6;k<9;k++)
00261       {
00262         if (nhit[k] != 0)
00263         {
00264           if (nhit[k] == 1)
00265           {
00266            super_hot[k] = 1;
00267           }
00268         }
00269 
00270         if (super_hot[k])
00271         {
00272           nsuper_hot_tel++;
00273         }
00274       }
00275 
00276 
00277       //If track, look in m2 prototype
00278       if (nsuper_hot_tel == 3)
00279       {
00280         //GetRMS(1) == 0 sigifie que tous les points ont la mewme valeur sur l'axe 1
00281         align = (hxy_tel_pad_temp->GetRMS(1) == 0 && hxy_tel_pad_temp->GetRMS(2) == 0); // GetRMS(1) RMS sur l'axe 1 ( x ) et GetRMS(2)  RMS sur l'axe 2 ( y )
00282 
00283         if (align)
00284         {
00285           if ( clusterInStrip[X] && nEntries[X] < MAXSTRIPENTRY )
00286           {
00287             strip[X][nEntries[X]] = (int)xStripNum;
00288             pad[X][nEntries[X]] = (Int_t)(pad_x/PADWIDTH);
00289             nEntries[X]++;
00290           }
00291 
00292           if ( clusterInStrip[Y] && nEntries[Y] < MAXSTRIPENTRY )
00293           {
00294             strip[Y][nEntries[Y]] = (int)yStripNum;
00295             pad[Y][nEntries[Y]] = (Int_t)(pad_y/PADWIDTH);
00296             nEntries[Y]++;
00297           }
00298           //Extrapolate track to m2 prototype
00299 
00300           for (int j=0;j<nchannel;j++)
00301           {
00302 
00303             MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00304             chbid = channel->GetChamberId();
00305 
00306             if (chbid == 10)
00307             {
00308 
00309               t2 = channel->GetBcIdDif();
00310               t3 = channel->GetBcIdHit();
00311               dt = t2 - t3;
00312 
00313               xpos = channel->GetX();
00314               ypos = channel->GetY();
00315               chipid = channel->GetChipId();
00316 
00317               if (dt>=dt_min && dt<=dt_max)
00318               {
00319                 hx_cor->Fill(pad_x,xpos);
00320                 hy_cor->Fill(pad_y,ypos);
00321 
00322                 hx_res->Fill(xpos - pad_x);
00323                 hy_res->Fill(ypos - pad_y);
00324 
00325               } // end if (dt>=dt_min && dt<=dt_max)
00326             } // end if (chbid == 10)
00327 
00328           } // end (int j=0;j<nchannel;j++)
00329 
00330         } // end if align
00331       } // end if (nsuper_hot_tel == 3)
00332 
00333       if ( i != nEvent-2 )
00334       {
00335          hxy_tel_pad_temp->Reset(); 
00336       }
00337     } // end for (int i=0;i<nEvent-1;i++)
00338   }  // end while TKey
00339 
00340 
00341   TCanvas c1;
00342   c1.Divide(2,3);
00343   c1.cd(1);
00344   cout << " n entries X "  << nEntries[X] << endl;
00345   TGraph *gx = new TGraph(nEntries[X],strip[X],pad[X]);
00346   gx->SetTitle ( "Correlation pad - strip X" ); 
00347   gx->Draw("A*");
00348   gx->Fit("pol1");
00349  
00350   cout << " n entries Y "  << nEntries[Y] << endl;
00351   TGraph *gy = new TGraph(nEntries[Y],strip[Y],pad[Y]);
00352   gy->SetTitle ( "Correlation pad - strip Y" ); 
00353   c1.cd(2);
00354   gy->Draw("A*");
00355   gy->Fit("pol1");
00356 
00357   c1.cd(3);
00358   hx_cor->Draw();
00359   hx_cor->ProfileX()->Fit("pol1");
00360   
00361   c1.cd(4);
00362   hy_cor->Draw();;
00363   hy_cor->ProfileX()->Fit("pol1");
00364 
00365   c1.cd(5);
00366   hx_res->Draw();
00367 
00368   c1.cd(6);
00369   hy_res->Draw();
00370 
00371 
00372   theApp->Run();
00373   delete theApp;
00374 
00375   return 0;
00376 }

Defines
#define	STRIPWIDTH 250
#define	PADWIDTH 10000
#define	MAXSTRIPENTRY 10000
#define	X 0
#define	Y 1
Functions
int	main (int argc, char **argv)
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/geomVerif.cpp File Reference

Defines

Functions

Define Documentation

Function Documentation
