Micromegas LAPP reconstruction

int main	(	int	argc,
		char **	argv
	)
Definition at line 40 of file yannis.cpp.
00040                               {
00041 
00042    if ( argc !=2  ) {
00043    FILE_LOG(logERROR)  << "usage: example rootFile " << endl;
00044    exit(1);
00045   }
00046 
00047    int search = 1;
00048    TH2F * hxy = new TH2F("hxy","",100,0,100,100,0,100);
00049    TH1I * hnhit = new TH1I("hnhit","Number of all hits",200,0,1000);
00050    TH1F * hDensity = new TH1F("hDensity","Hits Density ",200,0,50.);
00051    TH1I * hnhit_roi = new TH1I("hnhit_roi","Number of hits in tight DeltaT",130,0,130);
00052    TH1I * hnHitForw = new TH1I("hnhitForw","MIP Number of hits Forw",50,0,50);
00053    TH2F * hxy33 = new TH2F("hxy33","Chambre 34 Map",100,0,100,100,0,100);
00054    TH2F * hxy41 = new TH2F("hxy41","Chambre 49 Map",100,0,100,100,0,100);
00055    TH2F * hxy49 = new TH2F("hxy49","Chambre 49 Map",100,0,100,100,0,100);
00056    TProfile * tp_nhit = new TProfile("tp_nhit","",60,0,60,0,200,"");
00057    TH1F * hChi2X = new TH1F("hChi2X"," Chi2 X ", 200, 0., 10000.);
00058    TH1F * hChi2Y = new TH1F("hChi2Y"," Chi2 Y ", 200, 0., 10000.);
00059    TH1F * hResidY = new TH1F("hResidY"," Residual Y ", 200, -100., 100.);
00060    TH1F * hResidX = new TH1F("hResidX"," Residual X ", 200, -100., 100.);
00061    TProfile * pResidX = new TProfile("pResidX","Average Residual per layer X",50,0,50,-1000.,1000,"");
00062    TProfile * pResidY = new TProfile("pResidY","Average Residual per layer Y",50,0,50,-1000.,1000,"");
00063    TH1F * hPenteY = new TH1F("hPentY"," Slope Y ", 200, -0.1, 0.1);
00064    TH1F * hPenteX = new TH1F("hPentX"," Slope X ", 200, -0.1, 0.1);
00065    TH2F * hxyLayer0 = new TH2F("hxyLayer0","X verus Y fitted track. Layer 0",100,0,1000,100,0,1000);
00066 
00067    TH1I * hz = new TH1I("hz","",100,0,100);
00068 
00069    string rootName;
00070    rootName.assign(argv[1]);
00071    int nbHit = 0;
00072    int nbEvt = 0;
00073 
00074    int NbHitLayer[100];
00075    bool IsPenetratingMIP;
00076 
00077    float Zlayer[maxHits];
00078    float Xlayer[maxHits];
00079    float Ylayer[maxHits];
00080    float ErPos[maxHits];
00081   
00082   TFile f(rootName.c_str());
00083   TIter nextkey(f.GetListOfKeys());
00084   TKey *key;
00085 
00086  TFile tf("debug_Pions8GeV.root","RECREATE");
00087  int NbTrees = 0;
00088   while (key = (TKey*)nextkey()) 
00089   {
00090 
00091     TString keyname = key->GetName();
00092     if ((keyname == "hResidualX") || (keyname == "hResidualY")){}
00093     else{
00094 
00095     TTree *tree = (TTree*)key->ReadObj();                
00096     cout << "---------------------------- NEW TTREE -----------------"<< endl;
00097     nbEvt =0;
00098     nbHit = 0;
00099     NbTrees++;
00100     if( NbTrees > maxTrees ) continue;
00101     CaloRun *run=NULL;
00102 
00103     try {
00104       CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First());
00105       if ( run !- NULL )
00106       {
00107         cout << "Run[" << run->GetRunId() <<  "], temperature[" << run->GetTemperature() <<  "] pression[" << run->GetPressure() << "] " << endl ; 
00108       }
00109     }
00110     catch ( ... ) {}
00111     int RunId=1;
00112     TString name = Form("deltaTRun_%i",RunId);
00113     TString title = Form("delta distribution Run %i",RunId);
00114     TH1I * deltaTRun = new TH1I(name,title,20,-10,10);
00115 
00116     CaloEvent *evt =  new CaloEvent();
00117     TBranch *branch= tree->GetBranch("CaloEvent");
00118     branch->SetAddress(&evt);
00119 
00120     CaloHit* hit =NULL;
00121     int nbEntries = tree->GetEntries();
00122     
00123     for ( int evtNum = 0; evtNum < nbEntries ; evtNum++)
00124     {
00125       tree->GetEntry(evtNum);
00126       nbEvt++;
00127       if(nbEvt>maxEvents) continue;
00128       /*
00129       cout << "**** Event Information "<< dec << evtNum+1 << "/" << nbEntries<< " ****" <<  endl;
00130       cout << "    event num in TTree:" << evtNum <<endl;
00131       cout << "    event id:" << evt->fEventId <<endl;
00132 
00133       cout << "    number of channel hit:" << evt->GetNHits() <<endl;
00134       */
00135       int nbHits = evt->GetNHits();
00136       //      cout << "nb hits:" << nbHits << endl;
00137 
00138       hnhit->Fill(nbHits);
00139 
00140 
00141       if (nbHits > 120) continue;
00142 
00143 
00144       //number of hits in region of interest
00145       int nhit_roi = 0;
00146       int nhit_fit=0;
00147       float xpos = 0;
00148       float ypos = 0;
00149       float zpos=0;
00150       int deltat=0;
00151 
00152       for (int layer=0 ; layer <100 ; layer++){
00153         NbHitLayer[layer]=0;}
00154       for (int i=0 ; i<maxHits ;i++){
00155         Zlayer[i]=0;
00156         Xlayer[i]=0;
00157         Ylayer[i]=0;
00158         ErPos[i] = 30./sqrt(12.);
00159       }
00160 
00161       //Look for the Max Deltat
00162       //Find the 1st deltat and Book accordingly a histogram in order to find the maximum deltaT
00163         hit = (CaloHit*)evt->GetHits()->UncheckedAt(1);
00164         int deltat1 = hit->GetDeltaT();
00165         TH1I * hEvtDeltat = new TH1I("hEvtDeltat","Event DeltaT",20,deltat1-10,deltat1+10);
00166 
00167       for (int i=0 ; i <nbHits ; i++){
00168         hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);
00169         deltat=hit->GetDeltaT();
00170         hEvtDeltat->Fill(deltat);
00171       }
00172       int maxDeltat = deltat1-10 + hEvtDeltat->GetMaximumBin()-1;
00173       //      cout << "Max DeltaT " << maxDeltat << endl;
00174       //     hEvtDeltat->Write();
00175             delete hEvtDeltat;
00176 
00177       for(int i=0 ; i < nbHits ; i++)
00178       {
00179         hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);
00180 
00181         deltat=hit->GetDeltaT();
00182         deltaTRun->Fill(maxDeltat-deltat);
00183         if(abs(deltat-maxDeltat)>2) continue;
00184 
00185         xpos = hit->GetX()/10000+48;
00186         ypos = hit->GetY()/10000+48;
00187 
00188         zpos = hit->GetZ();
00189         if (zpos == 1370000){zpos = 1372000;}
00190         zpos = zpos / 28000;
00191 
00192         hz->Fill(zpos);
00193 
00194         //cout  << zpos << " " << " Y " << ypos <<endl;
00195 
00196         NbHitLayer[(int)zpos]++; 
00197 
00198         hxy->Fill(xpos,ypos);
00199         if((int)zpos == 33 ){hxy33->Fill(xpos,ypos);}
00200         if((int)zpos == 41 ){hxy41->Fill(xpos,ypos);}
00201         if((int)zpos == 49 ){hxy49->Fill(xpos,ypos);}
00202         //      cout << "deltat "<< deltat << endl;       
00203 
00204         nhit_roi++;
00205         // Prepare for fitting a muon
00206         // Choose layers to make the fit
00207         if (nhit_fit < maxHits){
00208           if(NbHitLayer[(int)zpos] >1 ) continue ; // keep only 1st hit in layers
00209           Zlayer[nhit_fit] = zpos*28; // Convert coordinates to mm
00210           Xlayer[nhit_fit] = xpos*10.;
00211           Ylayer[nhit_fit] = ypos*10.;
00212           nhit_fit++;
00213         }
00214           
00215       } // end hit
00216       hnhit_roi->Fill(nhit_roi);
00217 
00218       // Compute Hit Density to find MIP events
00219       int NbLayers=0;
00220       for (int layer=0 ; layer < 50 ; layer++){
00221         if(NbHitLayer[layer]>0)
00222           {
00223             NbLayers++;
00224             tp_nhit->Fill(layer+1,NbHitLayer[layer]);     
00225           }
00226       }
00227 
00228       float HitDensity = 1. * nhit_roi/NbLayers;
00229 
00230       //          cout << "Event num " << nbEvt << " Nb of Hits " << nbHits << " nb Layers " << NbLayers <<  " Density " << HitDensity << endl; 
00231 
00232       hDensity->Fill(HitDensity);
00233 
00234       //Look for penetrating MIPS
00235 
00236       int NbLayersForw=0;
00237       int NbLayersBack=0;
00238 
00239       if(1 < HitDensity < 3){
00240 
00241         for (int layer=0 ; layer < 10 ; layer++){
00242           //cout<<layer<<" "<<NbHitLayer[layer]<<endl;
00243           if(0 < NbHitLayer[layer] && NbHitLayer[layer] <= 3){NbLayersForw++;}
00244         }
00245         for (int layer=35 ; layer < 45 ; layer++){
00246           //cout<<layer<<" "<<NbHitLayer[layer]<<endl;
00247           if(0 < NbHitLayer[layer]  && NbHitLayer[layer] <= 3){NbLayersBack++;}
00248         }
00249 
00250         //cout<<NbLayersForw<<" "<<NbLayersBack<<endl;
00251 
00252         hnHitForw->Fill(NbLayersForw);
00253         hnHitForw->Fill(NbLayersBack+20);
00254 
00255         //Define penetrating MIP
00256         IsPenetratingMIP = 0;
00257         if(NbLayersForw>6 && NbLayersBack > 6){ IsPenetratingMIP = 1;}
00258       } // end Penetrating MIPS
00259 
00260 
00261       // Fit a straight line for penetrating MIPs 
00262        if(IsPenetratingMIP){
00263          //      for (int i=0;i<nhit_fit;i++){cout << i << " Z " << Zlayer[i] << " Y "<< Ylayer[i] << " X " << Xlayer[i] << endl;}
00264          float zmax = 50.* 28.;
00265         TGraphErrors *grx = new TGraphErrors(nhit_fit,Zlayer,Xlayer,0,ErPos);
00266         TF1 *f1 = new TF1("f1","pol1",0,zmax);
00267         grx->Fit("f1","Q");
00268         float betaX = f1->GetParameter(0);
00269         float alfaX = f1->GetParameter(1);
00270         double chi2X = f1->GetChisquare();
00271         TGraphErrors *gry = new TGraphErrors(nhit_fit,Zlayer,Ylayer,0,ErPos);
00272         TF1 *f2 = new TF1("f2","pol1",0,zmax);
00273         gry->Fit("f2","Q");
00274         float betaY = f2->GetParameter(0);
00275         float alfaY = f2->GetParameter(1);
00276         double chi2Y = f2->GetChisquare();
00277         hChi2X->Fill(chi2X);
00278         hChi2Y->Fill(chi2Y);
00279         hxyLayer0->Fill(betaX,betaY); // Make X Y Map at the first layer
00280 
00281         //              cout << " Pente X" << alfaX << "  offset " << betaX << " chi2 " << chi2X << endl;
00282         //              cout << " Pente Y" << alfaY << "  offset " << betaY << " chi2 " << chi2Y << endl;
00283 
00284         // Compute residuals to fitted track
00285         for (int j=1; j<nhit_fit ; j++){
00286           float ResX = Xlayer[j]-( alfaX*Zlayer[j]+betaX );
00287           float ResY = Ylayer[j]-( alfaY*Zlayer[j]+betaY );
00288           hResidY->Fill(ResY);
00289           hResidX->Fill(ResX);
00290           pResidX->Fill(Zlayer[j]/28,ResX);
00291           pResidY->Fill(Zlayer[j]/28,ResY);
00292         }
00293       }
00294  
00295       if (nbEvt%1000 == 0) {
00296       cout << "Processing event: " << nbEvt
00297            << endl;  
00298       }  
00299     } // end event
00300 
00301     deltaTRun->Write();
00302 
00303     
00304 
00305     } // end run ( TTree )
00306   
00307     
00308 }
00309   tf.cd();
00310   tp_nhit->Write();
00311   hxy->Write();
00312   hxy33->Write();
00313   hxy41->Write();
00314   hxy49->Write();
00315   hz->Write();
00316   hnhit->Write();
00317   hnhit_roi->Write();
00318   hDensity->Write();
00319   hnHitForw->Write();
00320   hChi2X->Write();
00321   hChi2Y->Write();
00322   hResidY->Write();
00323   hResidX->Write();
00324   pResidX->Write();
00325   pResidY->Write();
00326   hxyLayer0->Write();
00327 }

Functions
int	main (int argc, char **argv)
Variables
const int	maxEvents = 10000
const int	maxTrees = 2
const int	maxHits = 200
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/yannis.cpp File Reference

Functions

Variables

Function Documentation

Variable Documentation
