StRoot: StEEmcPool/StEzSmdCal/EEsmdCal.cxx Source File

00001 // $Id: EEsmdCal.cxx,v 1.20 2009/02/04 20:33:22 ogrebeny Exp $
00002  
00003 #include <assert.h>
00004 #include <stdlib.h>
00005 #include <string.h>
00006 
00007 #include <TClonesArray.h>
00008 #include <TObjArray.h> 
00009 #include <TH1.h> 
00010 #include <TH2.h> 
00011 #include <TFile.h> 
00012 
00013 #include "EEsmdCal.h"
00014 #include "EEsmdPlain.h"
00015 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00016 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
00017 #include "StEEmcUtil/EEmcSmdMap/EEmcSmdMap.h"
00018 #include "StEEmcUtil/database/cstructs/eemcConstDB.hh"
00019 #include "StEEmcUtil/database/EEmcDbItem.h"
00020 
00021 #ifdef StRootFREE
00022   #include "EEmcDb/EEmcDb.h"
00023 #else
00024   #include "StEEmcUtil/database/StEEmcDb.h"
00025 #endif
00026 
00027  
00028 ClassImp(EEsmdCal)
00029 //--------------------------------------------------
00030 //--------------------------------------------------
00031 EEsmdCal::EEsmdCal(){
00032   nInpEve=0; 
00033   HList=0; 
00034   eeDb=0;
00035 
00036   // clear all histo pointers just in case
00037   memset(hT,0,sizeof(hT));
00038   memset(hSs,0,sizeof(hSs));
00039   // memset(hSp,0,sizeof(hSp));
00040 
00041   dbMapped=-1;
00042   memset(dbT,0,sizeof(dbT));
00043   memset(dbS,0,sizeof(dbS));
00044 
00045   // initialization
00046   smdHitPl=new EEsmdPlain [MaxSmdPlains];
00047   geoTw=new EEmcGeomSimple;
00048   geoSmd= EEmcSmdGeom::instance();
00049   mapSmd= EEmcSmdMap::instance();
00050 
00051   thrMipSmdE=-1; emptyStripCount=-2; 
00052   twMipRelEneLow=-3; twMipRelEneHigh=-4;
00053   offCenter=0.7; thrMipPresAdc=-5;
00054   
00055   maxStripAdc=120; // suppress large jump in ped or sticky bits
00056 
00057   // chose which stat bits are fatal
00058   killStat=EEMCSTAT_ONLPED  | EEMCSTAT_HOTSTR ;
00059 
00060   printf("EEsmdCal() constructed\n");
00061 }
00062 
00063 //--------------------------------------------------
00064 //--------------------------------------------------
00065 EEsmdCal::~EEsmdCal() {/* noop */}
00066 
00067 
00068 //-------------------------------------------------
00069 //-------------------------------------------------
00070 void EEsmdCal::init( ){
00071   printf("EEsmdCal() init , calibrate sector=%d\n",sectID);
00072 
00073   initAuxHisto();
00074 
00075   initTileHistoAdc('a',"inclusive ADC ",kBlack);
00076   initTileHistoAdc('d',"ADC, tag: UxV 2*thr",kBlack);
00077   initTileHistoAdc('e',"ADC, tag: mipT UxV 2*thr",kBlue);
00078 
00079   initTileHistoEne('f',"energy, tag: UxV 2*thr",kBlack);
00080   initTileHistoEne('g',"energy, tag: mipT UxV 2*thr",kBlue);
00081  
00082   initSmdHist('a',"inclusive ADC");
00083   initSmdHist('b',"ADC, tag: best MIP",kBlue);
00084   initSmdEneHist('e',"Energy (K)+(K+1)*tgh(eta): best MIP",kBlack);
00085 
00086   //.................... initialize MIP finding algo for SMD
00087   int i;
00088   for(i=0;i<MaxSmdPlains;i++) {
00089     smdHitPl[i].set(thrMipSmdE,emptyStripCount,i+'U');
00090   }
00091 
00092   printf("use thrMipSmdE/MeV=%.2f emptyStripCount=%d  twMipRelEne/high=%.2f/%.2f offCenter=%.2f maxStripAdc=%.1f thrMipPresAdc=%d\n", thrMipSmdE*1000.,emptyStripCount,twMipRelEneLow, twMipRelEneHigh,offCenter,maxStripAdc,thrMipPresAdc);
00093   assert(sectID>0 && sectID<=MaxSectors);
00094 
00095   //....................... initilize MIP energy in towers
00096   assert(twMipRelEneLow< twMipRelEneHigh);
00097 
00098   for(int i=0;i<MaxEtaBins;i++){
00099     float etaValue=(geoTw->getEtaMean(i));
00100     float tghEta=TMath::TanH(etaValue);
00101     twTghEta[i]=tghEta;
00102     towerMipE[i]= 1./(2.89*tghEta); // GeV EM
00103     presMipE[i]=0.0009/tghEta; //GeV MIP, 5mm*1.8mm
00104   }
00105   
00106   smdAvrMipE=0.0013; //GeV MIP;  7mm*1.8 MeV/cm
00107   
00108 #if 0  //smdMap - histos for finding mapping
00109   int ij;
00110   for( ij=0;ij<12;ij++) {
00111     char t1[100];
00112     sprintf(t1,"hM%c",'a'+ij);
00113     hM[ij]=new TH1F(t1,t1,600,0.5,600.5);
00114     HList->Add( hM[ij]);
00115  }
00116 #endif
00117 
00118 }
00119 
00120 
00121 //-------------------------------------------------
00122 //-------------------------------------------------
00123 void EEsmdCal::initRun(int runID){
00124   printf(" EEsmdCal::initRun(%d)\n",runID);
00125 
00126 #if 1   //smdMap verification
00127   if(dbMapped>0)  {
00128     printf(" EEsmdCal::initRun(%d) N-th time, Ignore\n",runID);
00129     return; 
00130   }
00131 #endif
00132 
00133   assert(dbMapped<0); // at the moment DB reloading is not implemented/tested,JB
00134   mapTileDb();
00135   dbMapped=runID;
00136 
00137   // now init all what relays on DB inofrmation
00138   addTwMipEbarsToHisto(kGreen,'g');
00139 
00140   addPresMipEbarsToHisto(kGreen,'P');
00141   addPresMipEbarsToHisto(kGreen,'Q');
00142   addPresMipEbarsToHisto(kGreen,'R');
00143 
00144   addSmdMipEbarsToHisto(kGreen,'U');
00145   addSmdMipEbarsToHisto(kGreen,'V');
00146   histoGains();
00147 }
00148 
00149 //--------------------------------------------------
00150 //--------------------------------------------------
00151 void EEsmdCal:: mapTileDb(){
00152   printf("EEsmdCal:: mapTileDb()\n");
00153 
00154   //....................  Tower like ....................
00155   char cT[mxTile]={'T','P','Q','R'};
00156   int iT=0;
00157   for(iT=0;iT<mxTile;iT++) {
00158     for(char iSub=0; iSub<MaxSubSec; iSub++){
00159       for(int iEta=0; iEta<MaxEtaBins; iEta++){
00160         int iPhi=iSect*MaxSubSec+iSub;
00161         dbT[iT][iEta][iPhi]=eeDb->getTile(sectID,iSub+'A',iEta+1,cT[iT]);
00162         //printf("%d %d %d '%s' \n",iT,iEta,iPhi,dbT[iT][iEta][iPhi]->name);
00163       }
00164     }
00165   }
00166 
00167   //................ SMD strips
00168   int iu;
00169   for(iu=0;iu<MaxSmdPlains;iu++) {
00170     int istr;
00171     for(istr=0;istr<MaxSmdStrips;istr++) {
00172       if(istr>=288) break; // ugly, but I gave up, JB
00173       dbS[iu][istr]=eeDb->getByStrip(sectID,iu+'U',istr+1);
00174     }
00175   }
00176 
00177 }
00178 
00179 
00180 //-------------------------------------------------
00181 //-------------------------------------------------
00182 void EEsmdCal::clear(){ // called for every event
00183     memset(tileAdc,0,sizeof(tileAdc));
00184     memset(tileEne,0,sizeof(tileEne));
00185     memset(tileThr,0,sizeof(tileThr));
00186     memset(smdAdc,0,sizeof(smdAdc));
00187     memset(smdEne,0,sizeof(smdEne));
00188 
00189     memset(killT,true,sizeof(killT));// default is dead
00190 
00191     int i;
00192     for(i=0;i<MaxSmdPlains;i++) {
00193     smdHitPl[i].clear();
00194    }
00195 
00196 }
00197 //-------------------------------------------------
00198 //-------------------------------------------------
00199 void EEsmdCal::findSectorMip( ){ 
00200   /*...................................
00201     ....                           ....
00202     ....   main physics analysis   ....
00203     ................................... */
00204 
00205   hA[9]->Fill(1);
00206   fillSmdHisto_a();  // inclusive SMD histos 
00207 
00208   for(char iSub=0; iSub<MaxSubSec; iSub++){
00209     for(int iEta=0; iEta<MaxEtaBins; iEta++){      
00210       int iPhi=iSect*MaxSubSec+iSub;
00211       fillOneTailHisto('a', iEta,iPhi);   // inclusive histos 
00212     }
00213   }
00214   
00215   // searching for MIP in SMD
00216   int ret=getUxVmip();
00217 
00218   if(ret>1)  hA[9]->Fill(2);// counts multiple MIP's per both planes
00219   if(ret==1) hA[9]->Fill(3);// counts multiple MIP's per both planes
00220   
00221   int nU,nV;
00222   EEsmdPlain *plU=smdHitPl+0;
00223   EEsmdPlain *plV=smdHitPl+1;
00224 
00225 
00226 #if 0  //smdMap verification
00227   // verify mapping for a subset of strips
00228   // U-plan must have MIP
00229   if(plU->nMatch>0) {
00230     int ist;
00231     for(ist=0;ist<12;ist++)
00232     scanSpike(smdAdc[1][30+ist-1],hM[ist]);
00233   }
00234 #endif
00235 
00236   // calibrate P,Q,R with UxV
00237   for(nU=0;nU<plU->nMatch;nU++){
00238     for(nV=0;nV<plV->nMatch;nV++){
00239       hA[9]->Fill(4);// any UxV pair
00240       calibAllwithMip(plU->iStrip[nU],plV->iStrip[nV]);
00241     }
00242   }
00243 
00244   return;
00245 }
00246 
00247 
00248 //-------------------------------------------------
00249 //-------------------------------------------------
00250 void EEsmdCal::calibAllwithMip(int iStrU, int iStrV){
00251 
00252   // find MIP for a UxV pair of strips
00253 
00254   //printf("\n jj iSect=%d,iStrU=%d,iStrV=%d\n",iSect,iStrU,iStrV);
00255   TVector3 r=geoSmd->getIntersection (iSect,iStrU,iStrV);
00256  
00257   //printf(" UxV = %f %f %f\n", r.x(),r.y(),r.z());
00258  
00259   int     iSecX, iSubX, iEtaX;
00260   Float_t dphi, deta;
00261   int ret=geoTw->getTower(r, iSecX, iSubX, iEtaX,dphi, deta);
00262   //printf("ret=%d, isecX=%d isubX=%d, ietaX=%d dphi=%f, deta=%f\n",ret,iSecX, iSubX, iEtaX, dphi, deta);
00263 
00264   if(ret==0 || iSecX!=iSect) return;  
00265   //................ UxV is in a tower boundary within selected sector
00266 
00267   hA[9]->Fill(5);
00268   hA[10]->Fill(iStrU+1);
00269   hA[11]->Fill(iStrV+1);
00270   
00271   // select central MIP in tower
00272   int inCenter=(fabs(dphi)<offCenter && fabs(deta)<offCenter);
00273   ((TH2F*) hA[21])->Fill( r.x(),r.y());
00274 
00275   if(!inCenter) return; 
00276   //...................... assure central UxV hit in a tower
00277 
00278   hA[9]->Fill(6);
00279   int iPhiX=iSect*MaxSubSec+iSubX;
00280   hA[24]->Fill(iPhiX+iEtaX*MaxPhiBins);
00281 
00282   //................ auxuiliary variables 
00283   // ped corrected adc
00284   float adcT=tileAdc[kT][iEtaX][iPhiX];
00285   float adcP=tileAdc[kP][iEtaX][iPhiX];
00286   float adcQ=tileAdc[kQ][iEtaX][iPhiX];
00287   float adcR=tileAdc[kR][iEtaX][iPhiX];
00288 
00289   // calibrated energy
00290   float eneT=tileEne[kT][iEtaX][iPhiX];  // GeV
00291   float eneP=tileEne[kP][iEtaX][iPhiX]*1000; // MeV
00292   float eneQ=tileEne[kQ][iEtaX][iPhiX]*1000; // MeV
00293   float eneR=tileEne[kR][iEtaX][iPhiX]*1000; // MeV
00294 
00295   // logical conditions:
00296   // to recover few tiles mark all faild towers as with ADC>thres
00297   bool thrP=false,  thrQ=false,  thrR=false, thrDum=false;
00298   bool *thr_p[mxTile]={&thrDum,&thrP,&thrQ,&thrR};
00299 
00300   int iT;
00301   for (iT=kP; iT<=kR;iT++) {// loop over pre/post
00302     bool thr=tileThr[iT][iEtaX][iPhiX] ;   // .... ADC must be above ped
00303     thr= thr || killT[iT][iEtaX][iPhiX];   //.... account for masked pixels
00304 #if 0 //....  after pre/post calibration is avaliable
00305     float preMipRelEneLow=0.6, preMipRelEneHigh=3.0; // tmp
00306     float ene2Mip=tileEne[iT][iEtaX][iPhiX]/presMipE[iEtaX];
00307     thr= thr &&   ene2Mip>preMipRelEneLow &&  ene2Mip<preMipRelEneHigh; 
00308 #endif
00309     (*thr_p[iT])=thr;    // .... record final answer
00310   }
00311 
00312   //printf("bb1 (adc)  P=%.1f Q=%.1f R=%.1f iphiX=%d\n",tileAdc[1][iEtaX][iPhiX],tileAdc[2][iEtaX][iPhiX],tileAdc[3][iEtaX][iPhiX],iPhiX);
00313   //printf("bb2 (bool) thr P=%d Q=%d R=%d\n",thrP,thrQ,thrR);
00314 
00315 
00316   //printf("bb3 (bool) thr P=%d Q=%d R=%d\n",thrP,thrQ,thrR);
00317 
00318 
00319   // check MIP upper/lower limits
00320   float RelTwEne=tileEne[kT][iEtaX][iPhiX]/towerMipE[iEtaX];
00321   bool mipT=  RelTwEne>twMipRelEneLow &&  RelTwEne<twMipRelEneHigh;
00322   mipT=  mipT || killT[kT][iEtaX][iPhiX]; // recover dead tower
00323   // printf("iphi=%d ieta=%d Tene=%f mipEne=%f mipT=%d adcT=%.1f\n",iPhiX,iEtaX,tileEne[kT][iEtaX][iPhiX],towerMipE[iEtaX],mipT,tileAdc[kT][iEtaX][iPhiX]);
00324 
00325 
00326   if(thrR) hA[9]->Fill(7);
00327 
00328   if(mipT )hA[9]->Fill(8);
00329 
00330   //...........  MIP  ADC  spectra .........................
00331   int iCut='d'-'a';
00332   if(         thrP && thrQ && thrR ) hT[iCut][kT][iEtaX][iPhiX]->Fill(adcT);
00333   if( mipT &&         thrQ && thrR ) hT[iCut][kP][iEtaX][iPhiX]->Fill(adcP);
00334   if( mipT && thrP         && thrR ) hT[iCut][kQ][iEtaX][iPhiX]->Fill(adcQ);
00335   if( mipT && thrP && thrQ         ) hT[iCut][kR][iEtaX][iPhiX]->Fill(adcR);
00336 
00337 
00338   iCut='e'-'a';
00339   if( mipT && thrP && thrQ && thrR ){
00340     hT[iCut][kT][iEtaX][iPhiX]->Fill(adcT);
00341     hT[iCut][kP][iEtaX][iPhiX]->Fill(adcP);
00342     hT[iCut][kQ][iEtaX][iPhiX]->Fill(adcQ);
00343     hT[iCut][kR][iEtaX][iPhiX]->Fill(adcR);
00344   }
00345   
00346   // .................  calibrate (energy) MIP spectra ...............
00347   iCut='f'-'a';
00348   if(         thrP && thrQ && thrR ) hT[iCut][kT][iEtaX][iPhiX]->Fill(eneT);
00349   if( mipT &&         thrQ && thrR ) hT[iCut][kP][iEtaX][iPhiX]->Fill(eneP);
00350   if( mipT && thrP         && thrR ) hT[iCut][kQ][iEtaX][iPhiX]->Fill(eneQ);
00351   if( mipT && thrP && thrQ         ) hT[iCut][kR][iEtaX][iPhiX]->Fill(eneR);
00352 
00353   iCut='g'-'a';
00354   if( mipT && thrP && thrQ && thrR ){
00355     hT[iCut][kT][iEtaX][iPhiX]->Fill(eneT);
00356     hT[iCut][kP][iEtaX][iPhiX]->Fill(eneP);
00357     hT[iCut][kQ][iEtaX][iPhiX]->Fill(eneQ);
00358     hT[iCut][kR][iEtaX][iPhiX]->Fill(eneR);
00359     ((TH2F*) hA[22])->Fill( r.x(),r.y());
00360   }
00361 
00362   //..................... calibration of SMD strips
00363   if( mipT && thrP && thrQ && thrR ) {
00364  
00365     int iuv,i2;
00366     const int mx=2; // # of strips from given plain
00367     int iStr[MaxSmdPlains];
00368     iStr[0]=iStrU; // working pointers
00369     iStr[1]=iStrV;
00370     float eUV=0;// sum from both plains
00371     for(iuv=0;iuv<MaxSmdPlains;iuv++) {
00372       float e12=0;
00373       for(i2=0;i2<mx;i2++) {
00374         int istrip=iStr[iuv]+i2;
00375         float adc=smdAdc[iuv][istrip];
00376         float ene=smdEne[iuv][istrip]*twTghEta[iEtaX]*1000; // MeV, at normal angle     
00377         if(adc<2) continue;// drop pedestal, tmp
00378         e12+=ene; // sum pairs
00379         // re-calibratiop  of strips
00380         hSs['b'-'a'][iuv][istrip]->Fill(adc);
00381       }// end of loop over one plain
00382       eUV+=e12;
00383       // SMD energy for pairs
00384       int istrip1=iStr[iuv];
00385       ((TH2F*) hA[20])->Fill(istrip1+1,e12);
00386       hSs['e'-'a'][iuv][istrip1]->Fill(e12);      
00387       hA[14+iuv]->Fill(istrip1+1); 
00388 
00389     }// end of loop over UV plains
00390     ((TH2F*) hA[23])->Fill(iEtaX+1,eUV); 
00391   }
00392   
00393 }
00394 
00395 #if 0   
00396         // SMD light attenuation for pair of strips
00397         StructEEmcStrip* bar=geoSmd->getStripPtr(istrip,iuv,iSect);
00398         TVector3 r1=bar->end1;
00399         TVector3 r2=bar->end2;
00400         float dx,dy;
00401         if(r1.Eta()<r2.Eta()) {
00402           dx=r1.x()-r.x(); 
00403           dy=r1.y()-r.y();
00404         } else {
00405           dx=r2.x()-r.x();
00406           dy=r2.y()-r.y();
00407         } 
00408         float dist=sqrt(dx*dx+dy*dy);
00409         int iG=istrip/stripGang;
00410         assert(iG>=0 && iG <MaxAt);
00411         hSc[iuv][iSubX][iG]->Fill(ene); 
00412         hSeta[iuv][iSubX][iG]->Fill(r.Eta());
00413         hSdist[iuv][iSubX][iG]->Fill(dist);
00414 #endif
00415 
00416 //-------------------------------------------------
00417 //-------------------------------------------------
00418 int EEsmdCal::getUxVmip(){
00419   //  printf("\n\n EEsmdCal::getUxVmip() eve=%d\n",nInpEve);
00420 
00421   int iuv;
00422   for(iuv=0;iuv<MaxSmdPlains;iuv++) {
00423     EEsmdPlain *pl=smdHitPl+iuv;
00424     pl->scanAdc(smdEne[iuv], thrMipSmdE);// use if SMD gains are  known
00425     //pl->scanAdc(smdAdc[iuv], 15); // use if SMD gains are NOT known
00426     pl->findMipPattern();
00427     hA[12+iuv]->Fill(pl->nMatch);
00428     // pl->print(1);
00429     // printf("iuv=%d  nM=%d\n",iuv, pl->nMatch);  
00430   }
00431   
00432   int ret=smdHitPl[0].nMatch*smdHitPl[1].nMatch;
00433 
00434   //  printf("UxV ret=%d\n",ret);
00435   return ret;
00436 }
00437 
00438 
00439 #include <TLine.h> // just for drawing
00440 #include <TText.h> // just for drawing
00441 #include <TCanvas.h> // just for drawing
00442 
00443 //-------------------------------------------------
00444 //-------------------------------------------------
00445 void EEsmdCal::finish(int k){
00446   printf("\n  EEsmdCal::finish(sec=%d) nInpEve=%d \n",sectID,nInpEve);
00447 
00448   if(k<=0) return;
00449   // some test drawing:
00450   printf("drawing ...\n");
00451   //  TFile* f=new TFile("outSec5b/mip05b-8zAB.hist.root");
00452 
00453   int iuv=1;
00454   TString tt="sec"; tt+=sectID;
00455   TCanvas *c=new TCanvas(tt,tt,400,400);
00456   c->Divide(1,2);
00457   //  TH2F * h=(TH2F *)f->Get("xy05ct");
00458   TH2F * h=(TH2F *)hA[22];
00459 
00460   for(iuv=0;iuv<2;iuv++) {
00461     c->cd(iuv+1);
00462     h->Draw("colz");
00463     char uv='U'+iuv;
00464     int i;
00465     for(i=0;i<13;i++) {
00466       int strip=i*20+21;   
00467       StructEEmcStrip* bar=geoSmd->getStripPtr(strip-1,iuv,iSect);
00468       // printf("len=%f\n",bar->lenght);
00469       
00470       TVector3 r1=bar->end1;
00471       TVector3 r2=bar->end2;
00472       
00473       TLine *ln=new TLine(r1.x(),r1.y(),r2.x(),r2.y());
00474       ln->Draw();
00475       ln->SetLineColor(kRed);
00476 
00477       if(strip==161 ||strip==201 ||strip==241 ) continue;
00478       TString strT=uv; strT+=strip;
00479       TText *tx;
00480       if(uv=='U') 
00481         tx=new TText(r2.x()+1,r2.y()-5,strT);
00482       else {
00483         if(strip>140)
00484           tx=new TText(r2.x()-10,r2.y()-8,strT);
00485         else
00486           tx=new TText(r2.x()-20,r2.y(),strT);
00487       }
00488       tx->Draw();
00489     }
00490     char sub;
00491     for(sub='A'; sub<='E';sub++) {
00492       int i=sub-'A';
00493       float x=45 -i*8;
00494       float y=-55 -i*5;
00495       TString strT=sub;
00496       TText *tx=new TText(x,y,strT);
00497       tx->Draw();
00498     }
00499 
00500   }
00501   
00502 
00503 }
00504 
00505 
00506 #if 0 //smdMap verification
00507 //-------------------------------------------------
00508 //-------------------------------------------------
00509 void EEsmdCal:: scanSpike(float adc1, TH1F *h){
00510   float adcTh=20;
00511   if(adc1<adcTh) return;
00512   int istrip;
00513   for(istrip=0+3;istrip<MaxSmdStrips-3;istrip++) {
00514     int iuv;
00515     for(iuv=0;iuv<2;iuv++) {
00516 #if 0
00517       if(smdAdc[iuv][istrip-3]>=adcTh) continue;
00518       if(smdAdc[iuv][istrip-2]>=adcTh) continue;
00519       if(smdAdc[iuv][istrip+2]>=adcTh) continue;
00520       if(smdAdc[iuv][istrip+3]>=adcTh) continue;
00521 
00522       if(smdAdc[iuv][istrip-1]>=adcTh) continue;
00523       if(smdAdc[iuv][istrip+1]>=adcTh) continue;
00524 #endif
00525       if(smdAdc[iuv][istrip  ]< adcTh) continue;
00526       h->Fill(1+300*iuv+istrip);
00527     }
00528   }
00529 }
00530 #endif
00531