Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members  

EDAnalysis.cxx

Go to the documentation of this file.
00001 
00002 #include "TSystem.h"
00003 #include "TStyle.h"
00004 
00005 #include "bfemDisplay/EDAnalysis.h"
00006 #include "bfemDisplay/Enum.h"
00007 #include "bfemDisplay/EDHistos.h"
00008 
00009 
00010 TH1 *hACD[13];
00011 TH1 *hXGT[4];
00012 TH1 *hTKR[26];
00013 TH1 *hTKRtotL[26];
00014 TH1 *hTKRtotR[26];
00015 
00016 TH1F *LOGID;
00017 TH1F *NLOGS;
00018 TH1F *BESTRANGE;
00019 TH1F *ELOG;
00020 TH1F *ETOTADCCAL;
00021 TH1F *ETOTCAL;
00022 TH1F *ADC0;
00023 TH1F *ADC1;
00024 TH1F *ADC2;
00025 TH1F *ADC3;
00026 TH2F *ETOTVSNLOGS;
00027 
00028 
00029 EDHistos *h;
00030 
00031 
00032 
00033 
00034 EDAnalysis::EDAnalysis( EDControl *evtCtrl )
00035 
00036 {
00037 
00038   m_evtCtrl = evtCtrl;
00039 
00040   fTree = m_evtCtrl->getRawTree();
00041 
00042   event = m_evtCtrl->getEvent();
00043 
00044   h = m_evtCtrl->getHistos();
00045 
00046   m_ped = m_evtCtrl->getCALPeds();
00047 
00048   m_rawdata = m_evtCtrl->getRawDataType(); // Get type of input data: real or mc.
00049 
00050 
00051   m_StartEvent = 0;
00052 
00053 }
00054 
00055 
00056 void EDAnalysis::StartWithEvent(Int_t event) {
00057 
00058 // set a starting event for Go.
00059 
00060   m_StartEvent = event;
00061 
00062 }
00063 
00064 
00065 
00066 
00067 void EDAnalysis::Rewind()
00068 
00069 {
00070 
00071   //   Start input file at beginning
00072 
00073        m_StartEvent = 0;
00074 
00075        return;
00076 
00077 } 
00078 
00079 
00080 
00081 
00082 void EDAnalysis::Go(Int_t numEvents)
00083 
00084 {
00085   
00086   TH1F *hist;
00087   
00088   // Event loop_______________________________________ 
00089   
00090 
00091   fTree = m_evtCtrl->getRawTree();
00092 
00093   event = m_evtCtrl->getEvent();
00094 
00095   m_rawdata = m_evtCtrl->getRawDataType(); 
00096 
00097   if ( fTree == 0 ) return;
00098 
00099   if ( event ) event->Clean();
00100   
00101   
00102   Int_t nentries;
00103     
00104   nentries = m_evtCtrl->getEntries();
00105   
00106   
00107   printf("\nNum Events in File is: %i\n", nentries);
00108   
00109   Int_t curI;
00110   
00111   Int_t nMax = TMath::Min(numEvents+m_StartEvent,nentries);
00112   
00113   if (m_StartEvent == nentries) {
00114     
00115     printf(" all events in file read\n");
00116     
00117     return;
00118 
00119   }
00120 
00121 
00122 
00123   // Refresh histogram pointers_______________________________________
00124 
00125 
00126 
00127   TFile *histFile = (TFile*)gROOT->GetFile("Histograms.root");
00128 
00129 
00130   ADC0 = (TH1F*)histFile->Get("ADC0");
00131 
00132   ADC1 = (TH1F*)histFile->Get("ADC1");
00133 
00134   ADC2 = (TH1F*)histFile->Get("ADC2");
00135 
00136   ADC3 = (TH1F*)histFile->Get("ADC3");
00137 
00138   LOGID = (TH1F*)histFile->Get("LOGID");
00139 
00140   BESTRANGE = (TH1F*)histFile->Get("BESTRANGE");
00141 
00142   ELOG = (TH1F*)histFile->Get("ELOG");
00143 
00144   ETOTVSNLOGS = (TH2F*)histFile->Get("ETOTVSNLOGS");
00145 
00146   NLOGS = (TH1F*)histFile->Get("NLOGS");
00147 
00148   ETOTCAL = (TH1F*)histFile->Get("ETOTCAL");
00149 
00150   ETOTADCCAL = (TH1F*)histFile->Get("ETOTADCCAL");
00151 
00152   for (Int_t iacd = 0; iacd < NUM_ACD_HISTS; ++iacd) {
00153 
00154     hACD[iacd] = (TH1*)h->m_histList->At(iacd);
00155 
00156   }
00157 
00158   for (Int_t ixgt = 0; ixgt < NUM_XGT_HISTS; ++ixgt) {
00159 
00160     hXGT[ixgt] = (TH1*)h->m_histList->At(ixgt+13);
00161 
00162   }
00163 
00164   for (Int_t ilayer = 0; ilayer < NUM_TKR_HISTS; ++ilayer) {
00165 
00166     hTKR[ilayer] = (TH1*)h->m_histList->At(ilayer+28);
00167     
00168     hTKRtotL[ilayer] = (TH1*)h->m_histList->At(ilayer+54);
00169     
00170     hTKRtotR[ilayer] = (TH1*)h->m_histList->At(ilayer+80);
00171   }
00172 
00173    // Start analysis code _______________________________________
00174 
00175   m_Stop = kFALSE;
00176 
00177   for (Int_t ievent=m_StartEvent; ievent<nMax; ievent++, curI=ievent) {
00178 
00179 
00180     gSystem->ProcessEvents();
00181 
00182     if ( m_Stop ) break;
00183 
00184     if (event) event->Clean();
00185 
00186     if (ievent %100 == 0) printf(" Event: %i\n ", ievent);
00187 
00188     fTree->GetEvent(ievent);
00189 
00190    
00192     // ACD
00194 
00195     AcdTile *iTile;
00196 
00197     Int_t kADC = 1;
00198 
00199     if(kADC) {
00200 
00201       Int_t nACD = event->getAcdDigi()->GetEntries();
00202 
00203       Float_t eTotalACD;        
00204           
00205       for (Int_t ihit=0; ihit < nACD; ihit++) {   
00206  
00207         event->getAcdDigi()->Sort();
00208 
00209         AcdTile *hit = (AcdTile*)event->getAcdDigi()->At(ihit);
00210 
00211         AcdId *tile = (AcdId*) hit->getId();
00212 
00213         UInt_t id = tile->getId();
00214 
00215         if (id >1000) continue;
00216 
00217         iTile = ( AcdTile* )event->getAcdTile(id);
00218 
00219         eTotalACD = (Float_t) iTile->getPulseHeight();
00220 
00221         hist = (TH1F*)hACD[ihit];
00222 
00223         hist->Fill(eTotalACD);
00224 
00225       }
00226 
00227     }
00228 
00229 
00231     // XGT
00233 
00234 
00235     Int_t kXGT = 1;
00236 
00237     if(kXGT) {
00238 
00239       Int_t nXGT = event->getXgtDigi()->GetEntries();
00240           
00241       for (Int_t ihit=0; ihit < nXGT; ihit++) {   
00242  
00243         event->getXgtDigi()->Sort();
00244 
00245         AcdTile *hit = (AcdTile*)event->getXgtDigi()->At(ihit);
00246 
00247         AcdId *xgt = (AcdId*) hit->getId();
00248 
00249         Int_t id = xgt->getId();
00250 
00251         Float_t eTotalXGT;
00252 
00253         iTile = ( AcdTile* )event->getXgt(id);
00254 
00255         eTotalXGT = iTile->getPulseHeight();
00256 
00257         hist = (TH1F*)hXGT[ihit];
00258 
00259         hist->Fill(eTotalXGT);
00260 
00261       }
00262 
00263     }
00264 
00265 
00266 
00268 
00269     // Tracker...
00270 
00272 
00273 
00274 
00275     Int_t kTKR = 1;
00276 
00277     if(kTKR) {  
00278 
00279       event->getTkrDigi()->Sort();
00280 
00281       int nTKR = event->getTkrDigi()->GetEntries();
00282 
00283       for (Int_t ilay = 0; ilay < nTKR; ilay++) {
00284 
00285         const TkrLayer* layer = event->getTkrLayer(ilay); 
00286 
00287         // getXY = 0 measures Y ==> strips along X
00288         
00289         //      Int_t isY = (Int_t)layer->getXY();
00290 
00291         Int_t LayerNum= layer->getLayerNum();
00292        
00293  
00294         hist = (TH1F*)hTKR[LayerNum];       
00295             
00296         TObjArray  *stripList = 0;
00297             
00298         stripList = ( ( TkrLayer* )layer )->getStrips();
00299 
00300         Int_t nStrips = stripList->GetEntries();
00301 
00302         for (Int_t iStrip=0; iStrip < nStrips; iStrip++) {
00303 
00304           StripId* strip = (StripId*)(stripList->At(iStrip));
00305 
00306           Int_t stripID = strip->getId();
00307 
00308           if (stripID <=1600) hist->Fill(stripID);
00309 
00310         }
00311 
00312         Int_t totR =  ( ( TkrLayer* )layer )->getToT(0);
00313          
00314         Int_t totL =  ( ( TkrLayer* )layer )->getToT(1);
00315 
00316         hist = (TH1F*)hTKRtotL[LayerNum];
00317             
00318         if (totL) hist->Fill(totL);
00319 
00320         hist = (TH1F*)hTKRtotR[LayerNum];
00321             
00322         if (totR) hist->Fill(totR);
00323 
00324       }
00325       
00326     }
00327 
00328 
00329 
00331 
00332     // Calorimeter...
00333 
00335 
00336 
00337 
00338     Int_t kCAL = 1;
00339 
00340     if(kCAL) {  
00341 
00342       Int_t nCAL = event->getCalDigi()->GetEntries();
00343 
00344       Int_t nLOGS = 0;
00345 
00346       Float_t eLog = 0;
00347 
00348       Float_t eTotcal = 0, eTotadc = 0;
00349 
00350       Int_t index;
00351 
00352       //Float_t gain[4] = {0.000065, 0.00022, 0.0037, 0.022};
00353       Float_t gain[4] = {0.000075, 0.0003, 0.0043, 0.025};       //Eric's gain Sept. 11 2001  - in MeV
00354          
00355       
00356       for (Int_t ihit=0; ihit < nCAL; ihit++) {
00357 
00358         Int_t En = 1;
00359 
00360         CalLog *hitcal = (CalLog*) event->getCalDigi()->At(ihit);       
00361 
00362         LogId *log = (LogId*)hitcal->getLogId();
00363 
00364         Int_t id = log->getId(); // ==> id differs from the old id! range now from 0 to 39
00365 
00366         Int_t lay = log->getLayer();
00367 
00368         Int_t column = log->getColumn();
00369 
00370         Int_t iped = lay*10+column;
00371 
00372         Int_t hitADC;
00373 
00374 
00375         ((TH1*)h->m_histList->At(17))->Fill(id); // LOGID
00376            
00377            
00378 
00379         for (Int_t end=0; end <= 1; end++) {
00380              
00381                      
00382           Float_t eMax = -1000.; Int_t bestRange = -1;
00383 
00384           for (Int_t adc = 0; adc <= 3; adc++) {
00385 
00386             if (!end) hitADC = 4095 - hitcal->getAdc(hitcal->POS,adc);
00387             else hitADC = 4095 - hitcal->getAdc(hitcal->NEG,adc);
00388 
00389             if (hitADC < 1) continue;
00390 
00391             if(adc == 0) ((TH1*)h->m_histList->At(23))->Fill(hitADC); // ADC0
00392 
00393             if(adc == 1) ((TH1*)h->m_histList->At(24))->Fill(hitADC); // ADC1
00394 
00395             if(adc == 2) ((TH1*)h->m_histList->At(25))->Fill(hitADC); // ADC2
00396 
00397             if(adc == 3) ((TH1*)h->m_histList->At(26))->Fill(hitADC); // ADC3
00398 
00399             if (hitADC > 3900) continue;
00400 
00401             if (hitADC > eMax) {
00402 
00403               eMax = hitADC;
00404 
00405               bestRange = adc;
00406 
00407             }
00408 
00409           }
00410 
00411 
00412 
00413           if(!end) index = 2 * bestRange;
00414 
00415           else index = 2 * bestRange + 1;
00416 
00417           // Correct for unexpected ordering in ROOT ped file.
00418           
00419           /*      if(index==1) index=2;
00420           
00421                   if(index==2) index=1;*/
00422 
00423 
00424 
00425           ((TH1*)h->m_histList->At(19))->Fill(bestRange); // BESTRANGE
00426 
00427           Float_t eHit;
00428 
00429           if (bestRange > -1) {
00430 
00431             if ( m_rawdata ) {  // real raw data 
00432 
00433 
00434             if ( m_ped ) {
00435                
00436               TVector *vv = (TVector *)m_ped->At(index);
00437                                
00438               Float_t calped = vv->operator()( iped );
00439 
00440               eHit = (eMax-calped)*gain[bestRange];
00441 
00442               if (eHit < 0) En = 0;
00443 
00444             }
00445 
00446             else { eHit = 0; En = 0; }
00447             }
00448 
00449             else { // mc data
00450 
00451               eHit = eMax*gain[bestRange];
00452               if (eHit <= 0) En = 0;
00453 
00454             }
00455 
00456             eLog += eHit/2.;
00457 
00458             eTotadc += eMax/2.;       
00459 
00460 
00461           }
00462 
00463         }
00464 
00465         if (En) {
00466 
00467           ((TH1*)h->m_histList->At(20))->Fill(eLog); // ELOG
00468 
00469           eTotcal += eLog;
00470 
00471           nLOGS++;
00472 
00473         }
00474 
00475         eLog = 0;
00476 
00477       } // End CAL logs loop
00478 
00479          
00480 
00481       if (eTotcal > 0) ((TH1*)h->m_histList->At(22))->Fill(eTotcal); // ETOTCAL
00482       
00483       ((TH1*)h->m_histList->At(21))->Fill(eTotadc); // ETOTADCCAL
00484 
00485       ((TH1*)h->m_histList->At(18))->Fill(nLOGS); // NLOGS
00486 
00487       ((TH1*)h->m_histList->At(27))->Fill(eTotadc,(Float_t)nLOGS); // ETOTVSNLOGS
00488 
00489               
00490 
00491     }
00492 
00493 
00494 
00495   } // End event loop
00496 
00497   // end analysis code in event loop ___________________________
00498 
00499 
00500   histFile->Write();
00501 
00502   m_StartEvent = curI;
00503 
00504       
00505 
00506   cout << "\nAnalysis finished at event " << curI << ".\n";
00507 
00508 
00509 } // End GO()
00510 
00511 
00512 
00513 
00514 
00515 
00516 
00517 
00518 
00519 
00520 
00521 
00522 
00523 
00524 

Generated at Mon Nov 26 18:20:06 2001 by doxygen1.2.3 written by Dimitri van Heesch, © 1997-2000