StRoot: StEEmcPool/StEzSmdCal/EEsmdCalHisto.cxx Source File

00001 // $Id: EEsmdCalHisto.cxx,v 1.18 2009/12/03 22:35:03 ogrebeny Exp $
00002  
00003 #include <assert.h>
00004 #include <stdlib.h>
00005 
00006 #include <TClonesArray.h>
00007 #include <TObjArray.h> 
00008 #include <TH1.h> 
00009 #include <TH2.h> 
00010 #include <TLine.h> 
00011 #include <TFile.h> 
00012 
00013 #include "EEsmdCal.h"
00014 #include "StEEmcUtil/database/EEmcDbItem.h"
00015 
00016 #ifdef StRootFREE
00017   #include "EEmcDb/EEmcDb.h"
00018 #else
00019   #include "StEEmcUtil/database/StEEmcDb.h"
00020 #endif
00021 
00022 //--------------------------------------------------
00023 //--------------------------------------------------
00024 void EEsmdCal::initTileHistoAdc(char cut, const char *title, int col) {
00025   int iCut=cut -'a';
00026   assert(iCut>=0 && iCut<kCut);
00027   char tt1[100], tt2[500];
00028 
00029   const char *cTile[mxTile]={"Tower","Pres1","Pres2","Post"};
00030   const char cT[mxTile]={'T','P','Q','R'};
00031 
00032   int iT=0;
00033   for(iT=0;iT<mxTile;iT++) {
00034     for(char iSub=0; iSub<MaxSubSec; iSub++){
00035       for(int iEta=0; iEta<MaxEtaBins; iEta++){
00036         char sub=iSub+'A';
00037         int eta=iEta+1;
00038         int iPhi=iSect*MaxSubSec+iSub;
00039         char core[100];
00040         sprintf(core,"%02d%c%c%02d",sectID,cT[iT],sub,eta);
00041         sprintf(tt1,"%c%s",cut,core);
00042         sprintf(tt2,"%s(%c) %s , %s; ADC-ped",cTile[iT],cut,core,title);
00043         //printf("tt1=%s, tt2=%s\n",tt1,tt2);
00044         TH1F *h;
00045         if(cT[iT]=='T') 
00046           h=new TH1F(tt1,tt2,220,-20,200.);
00047         else 
00048           h=new TH1F(tt1,tt2,400,-50,350.);
00049         h->SetLineColor(col);
00050         HList->Add(h);
00051         hT[iCut][iT][iEta][iPhi]=h;
00052       }
00053     }
00054   }
00055 }
00056 
00057 
00058 
00059 
00060 //--------------------------------------------------
00061 //--------------------------------------------------
00062 void EEsmdCal::initTileHistoEne(char cut, const char *title, int col) {
00063   int iCut=cut -'a';
00064   assert(iCut>=0 && iCut<kCut);
00065   char tt1[100], tt2[500];
00066 
00067   const char *cTile[mxTile]={"Tower","Pres1","Pres2","Post"};
00068   const char cT[mxTile]={'T','P','Q','R'};
00069   const char *cUnits[mxTile]={"(GeV)_EM", "(MEV)_MIP", "(MEV)_MIP", "(MEV)_MIP"};
00070   
00071   int iT=0;
00072   for(iT=0;iT<mxTile;iT++) {
00073     for(char iSub=0; iSub<MaxSubSec; iSub++){
00074       for(int iEta=0; iEta<MaxEtaBins; iEta++){
00075         char sub=iSub+'A';
00076         int eta=iEta+1;
00077         int iPhi=iSect*MaxSubSec+iSub;
00078         char core[100];
00079         sprintf(core,"%02d%c%c%02d",sectID,cT[iT],sub,eta);
00080         sprintf(tt1,"%c%s",cut,core);
00081         sprintf(tt2,"%s(%c) %s , %s; ADC-ped/gain %s",cTile[iT],cut,core,title,cUnits[iT]);
00082         //printf("tt1=%s, tt2=%s\n",tt1,tt2);
00083         TH1F *h=new TH1F(tt1,tt2,200,-0.5,9.5);
00084         h->SetLineColor(col);
00085         HList->Add(h);
00086         hT[iCut][iT][iEta][iPhi]=h;
00087       }
00088     }
00089   }
00090 }
00091 
00092 
00093 
00094 //--------------------------------------------------
00095 //--------------------------------------------------
00096 void EEsmdCal::addTwMipEbarsToHisto (int col, char mxC) {
00097   // search all existing tower histo (with '05T' in name) and add 
00098   // bars for MIP limits
00099 
00100   assert(dbMapped>0);
00101 
00102   char core[100];
00103   sprintf(core,"%02dT",sectID);
00104   float yMax=1000;
00105   TIterator *it=HList->MakeIterator();
00106   TH1 *h;
00107   while( (h=(TH1*) it->Next())) {
00108     const char *name=h->GetName();
00109     if(strstr(name,core)==0) continue;
00110     if(name[0]>mxC) continue; // to skip energy-like histos declared later
00111     //printf("%s\n",h->GetName());
00112     int iSub=name[4]-'A';
00113     int iEta=atoi(name+5)-1;
00114     int iPhi=iSect*MaxSubSec+iSub;
00115     const EEmcDbItem *x=dbT[kT][iEta][iPhi];
00116     assert(x);
00117     if(x->gain<=0) continue;
00118     TList *L=h->GetListOfFunctions();
00119 
00120     float adcC=towerMipE[iEta]*x->gain;
00121     if(name[0]>'e') adcC=towerMipE[iEta]; //ugly hack,JB, last two not ADC
00122 
00123     TLine *ln=new TLine(adcC,0,adcC,yMax);
00124     ln->SetLineColor(kRed); ln->SetLineStyle(2);
00125     L->Add(ln);
00126 
00127     float adc=adcC*twMipRelEneLow;
00128     ln=new TLine(adc,0,adc,yMax);
00129     ln->SetLineColor(col);
00130     L->Add(ln);
00131 
00132     adc=adcC*twMipRelEneHigh;
00133     ln=new TLine(adc,0,adc,yMax);
00134     ln->SetLineColor(col);
00135     L->Add(ln);
00136   }
00137 }
00138 
00139 //--------------------------------------------------
00140 //--------------------------------------------------
00141 void EEsmdCal::addPresMipEbarsToHisto (int col, char cT) {
00142   // search all existing tower histo (with '05X' in name) and add 
00143   // bars for MIP limits
00144 
00145   int iT=1+cT-'P';
00146   assert(iT>kT && iT<mxTile);
00147   assert(dbMapped>0);
00148 
00149   char core[100];
00150   sprintf(core,"%02d%c",sectID,cT);
00151   float yMax=1000;
00152   TIterator *it=HList->MakeIterator();
00153   TH1 *h;
00154   while( (h=(TH1*) it->Next())) {
00155     const char *name=h->GetName();
00156     if(strstr(name,core)==0) continue;
00157     // printf("%s\n",h->GetName());
00158     int iSub=name[4]-'A';
00159     int iEta=atoi(name+5)-1;
00160     int iPhi=iSect*MaxSubSec+iSub;
00161     const EEmcDbItem *x=dbT[iT][iEta][iPhi];
00162     assert(x);
00163     if(x->gain<=0) continue;
00164     
00165     float adcC=presMipE[iEta]*x->gain;
00166     if(name[0]>'e') adcC=presMipE[iEta]*1000; //ugly hack2,JB, last two not ADC
00167 
00168     TLine *ln=new TLine(adcC,0,adcC,yMax);
00169     ln->SetLineColor(col); ln->SetLineStyle(2);
00170     TList *L=h->GetListOfFunctions();
00171     L->Add(ln);
00172   }
00173 }
00174 
00175 
00176 //--------------------------------------------------
00177 //--------------------------------------------------
00178 void EEsmdCal::initSmdHist(char cut, const char *title, int col) {
00179   int iCut=cut -'a';
00180   assert(iCut>=0 && iCut<kCut);
00181   char tt1[100], tt2[500];
00182 
00183   int iuv,istrip;
00184   for(iuv=0;iuv<MaxSmdPlains;iuv++) {
00185     for(istrip=0;istrip<MaxSmdStrips;istrip++) {
00186         char core[100];
00187         // single strip (ADC-ped) spectra
00188         sprintf(core,"%02d%c%03d",sectID,iuv+'U',istrip+1);
00189         sprintf(tt1,"%c%s",cut,core);
00190         sprintf(tt2,"SMD(%c) %s , %s; ADC-ped",cut,core,title);
00191         //printf("tt1=%s, tt2=%s\n",tt1,tt2);
00192         TH1F *h=new TH1F(tt1,tt2,400,-50,350);
00193         //TH1F *h=new TH1F(tt1,tt2,4400,-400,4000);//tmp
00194         h->SetLineColor(col);
00195 
00196         HList->Add(h);
00197         hSs[iCut][iuv][istrip]=h;
00198     }
00199   }
00200 
00201 }
00202 
00203 //--------------------------------------------------
00204 //--------------------------------------------------
00205 void EEsmdCal::initSmdEneHist(char cut, const char *title, int col) {
00206   int iCut=cut -'a';
00207   assert(iCut>=0 && iCut<kCut);
00208   char tt1[100], tt2[500];
00209 
00210   int iuv,istrip;
00211   for(iuv=0;iuv<MaxSmdPlains;iuv++) {
00212     for(istrip=0;istrip<MaxSmdStrips;istrip++) {
00213         char core[100]; 
00214         sprintf(core,"%02d%c%03d",sectID,iuv+'U',istrip+1);
00215         sprintf(tt1,"%c%s",cut,core);
00216         sprintf(tt2,"SMD(%c) %s , %s; normal incident MIP energy (MeV)",cut,core,title);
00217         //printf("tt1=%s, tt2=%s\n",tt1,tt2);
00218         TH1F *h=new TH1F(tt1,tt2,200,-0.1,9.9);
00219         h->SetLineColor(col);
00220 
00221         HList->Add(h);
00222         hSs[iCut][iuv][istrip]=h;
00223     }
00224   }
00225 
00226 }
00227 
00228 //--------------------------------------------------
00229 //--------------------------------------------------
00230 void EEsmdCal::addSmdMipEbarsToHisto (int col, char cU) {
00231   // search all existing tower histo (with 'x05U' in name) and add 
00232   // bars for MIP limits
00233   
00234   assert(cU=='U' || cU=='V');
00235   assert(dbMapped>0);
00236 
00237   char core[100];
00238   sprintf(core,"%02d%c",sectID,cU);
00239   float yMax=1000;
00240   TIterator *it=HList->MakeIterator();
00241   TH1 *h;
00242   while( (h=(TH1*) it->Next())) {
00243     const char *name=h->GetName();
00244     if(strstr(name,core)!=name+1) continue;
00245     // printf("%s-%c\n",h->GetName(),h->GetName()[0]);
00246     //printf("%s\n",h->GetTitle());
00247     int iU=name[3]-'U';
00248     int iStr=atoi(name+4)-1;
00249     assert(iU>=0 && iU<MaxSmdPlains);
00250     assert(iStr>=0 && iStr<MaxSmdStrips);
00251     const EEmcDbItem *x=dbS[iU][iStr];
00252     //   printf("iU=%c iStr=%d\n",iU,iStr);
00253     if(x==0) continue;
00254 
00255     if((h->GetName()[0])=='a') {// one more hack
00256       char tt3[500];
00257       sprintf(tt3,"%s, tube=%s",h->GetTitle(),x->tube);
00258       h->SetTitle(tt3);
00259       // printf("%s %p %s\n",h->GetTitle(),h,h->GetName());
00260     }
00261 
00262     if(x->gain<=0) continue;
00263     TList *L=h->GetListOfFunctions();    
00264     float adcC=smdAvrMipE*x->gain;
00265     if((h->GetName()[0])>'b')  adcC=smdAvrMipE*1000.; // now in MeV
00266 
00267     TLine *ln=new TLine(adcC,0,adcC,yMax);
00268     ln->SetLineColor(col); ln->SetLineStyle(2);
00269     L->Add(ln);
00270     
00271     // printf("%s iU=%c iStr=%d adcC=%f\n",name,iU+'U',iStr,adcC);
00272 
00273   }
00274 }
00275 
00276 
00277 //--------------------------------------------------
00278 //--------------------------------------------------
00279 void EEsmdCal::histoGains(){
00280 
00281   //.................... SMD ................
00282   int iuv,istrip;
00283   for(iuv=0;iuv<MaxSmdPlains;iuv++) {
00284     for(istrip=0;istrip<MaxSmdStrips;istrip++) {
00285       const EEmcDbItem *x=dbS[iuv][istrip];
00286       if(x==0) continue;
00287       //if(x->fail) continue; // use any non-zero gain
00288       if(x->gain<=0) continue;
00289       hA[16+iuv]->Fill(x->strip,x->gain);
00290       // dig out MAPMT pixel
00291       const char *tube=x->tube+2;
00292       assert(tube[0]=='S');
00293       int box=atoi(tube+1);
00294       int pmt=atoi(tube+3);
00295       int pix=atoi(tube+6);
00296       int xx=(pmt-1)*16+pix;
00297       hA[5+box-1]->Fill(xx,x->gain);
00298       //  printf("%d %d %d %d", box,pmt,pix,xx); x->print();      assert(1==3);
00299     }
00300   }
00301 
00302 
00303   //.............. Pre/post
00304   int iT,iEta,iPhi;
00305   for(iT=kP;iT<=kR;iT++) 
00306     for(iEta=0;iEta<MaxEtaBins;iEta++)
00307       for(iPhi=0;iPhi<MaxPhiBins;iPhi++){
00308         const EEmcDbItem *x=dbT[iT][iEta][iPhi];
00309         if(x==0) continue;
00310         // if(x->fail) continue;  // use any non-zero gain
00311         if(x->gain<=0) continue;
00312         // dig out the MAPMT pixel
00313         assert(x->sec==sectID);
00314         const char *tube=x->tube+2;
00315         assert(tube[0]=='P');
00316         int pmt=atoi(tube+3);
00317         int pix=atoi(tube+6);
00318         int xx=(pmt-1)*16+pix;
00319         hA[5+3]->Fill(xx,x->gain);
00320       }
00321 
00322 
00323   //.............. Towers
00324   iT=0;
00325   for(iEta=0;iEta<MaxEtaBins;iEta++)
00326     for(iPhi=0;iPhi<MaxPhiBins;iPhi++){
00327       const EEmcDbItem *x=dbT[iT][iEta][iPhi];
00328       if(x==0) continue;        
00329       if(x->gain<=0) continue;
00330       assert(x->sec==sectID);
00331       int ispir=(x->eta-1)*MaxSubSec +  x->sub-'A';
00332       hA[4]->Fill(ispir,x->gain);
00333     }
00334 
00335 }
00336 
00337 
00338 //--------------------------------------------------
00339 //--------------------------------------------------
00340 void EEsmdCal::initAuxHisto(){
00341   int i;
00342   //  float Emax=2;
00343   memset(hA,0,sizeof(hA));
00344 
00345   char tt1[100], tt2[500], tt0[100]; 
00346   TH1F *h;
00347   TH2F *h2;
00348 
00349   // tower gains
00350   sprintf(tt1,"ug%02dT",sectID);
00351   sprintf(tt2,"used tower gains sec=%d ; x=spiral=(eta-1)*5+subs-A; gain [ch/GeV]",sectID);
00352   h=new TH1F(tt1,tt2, 60,-0.5,59.5);
00353   hA[4]=h;    
00354 
00355   for(i=0;i<4;i++) {
00356     // gains vs. MAPMT box/pmt/ch, histo 5,6,7,8    
00357     if(i<3) { 
00358       sprintf(tt0,"%02d_S%d",sectID,i+1);
00359     } else {
00360       sprintf(tt0,"%02d_P1",sectID);
00361     }
00362     sprintf(tt1,"ug%s",tt0);
00363     sprintf(tt2,"used gains MAPMT box %s ; chann=(pmt ID-1)*16+pix ; gain [ch/GeV]",tt0);
00364 
00365     h=new TH1F(tt1,tt2, MaxMapmtCrateCh,0.5, MaxMapmtCrateCh+0.5);
00366     hA[5+i]=h;    
00367   }
00368 
00369   sprintf(tt1,"my%02dStat",sectID);
00370   hA[9]=new TH1F (tt1,"type of events ",30,.5,30.5);
00371 
00372   for(i=0;i<MaxSmdPlains;i++) {    
00373     sprintf(tt1,"fr%02d%c",sectID,i+'U');
00374     sprintf(tt2,"freq. of MIP, UxV only, plane %02d%c; strip ID",sectID,i+'U');
00375     h=new TH1F(tt1,tt2,MaxSmdStrips,-0.5,MaxSmdStrips-0.5);
00376     hA[10+i]=h;
00377     
00378     sprintf(tt1,"mm%02d%c",sectID,i+'U');
00379     sprintf(tt2,"freq of 00xx00  pattern per plane %02d%c",sectID,i+'U');
00380     h=new TH1F(tt1,tt2,20,-0.5,19.5);
00381     hA[12+i]=h;
00382 
00383     sprintf(tt1,"fr%02d%cm",sectID,i+'U');
00384     sprintf(tt2,"freq. of best MIP, plane %02d%c; strip ID",sectID,i+'U');
00385     h=new TH1F(tt1,tt2,MaxSmdStrips,-0.5,MaxSmdStrips-0.5);
00386     hA[14+i]=h;
00387     
00388     sprintf(tt1,"ug%02d%c",sectID,i+'U');
00389     sprintf(tt2,"used gains for plane %02d%c; strip ID; gain [ch/GeV]",sectID,i+'U');
00390     h=new TH1F(tt1,tt2,MaxSmdStrips,0.5,MaxSmdStrips+0.5);
00391     hA[16+i]=h;    
00392   }
00393 
00394   //..................
00395   sprintf(tt1,"ep%02dUorV",sectID);
00396   sprintf(tt2,"#Sigma E of 2-strips, normal angle, best MIP , SMD %02dUorV; strip ID; MIP #Sigma E (MeV) ",sectID);
00397   h2=new TH2F(tt1,tt2,30,0,300,100,-.1,10.);
00398   hA[20]=(TH1F*)h2;
00399   
00400   //..................
00401   sprintf(tt1,"xy%02d",sectID);
00402   sprintf(tt2,"MIP position , UxV only, sect=%02d; X(cm); Y(cm) ",sectID);
00403   h2=new TH2F(tt1,tt2,500,-250,250,500,-250,250);
00404   hA[21]=(TH1F*)h2;
00405   
00406   sprintf(tt1,"xy%02dm",sectID);
00407   sprintf(tt2,"MIP position, best MIP, sect=%02d; X(cm); Y(cm) ",sectID);
00408   h2=new TH2F(tt1,tt2,250,-250,250,250,-250,250);
00409   hA[22]=(TH1F*)h2; 
00410    
00411   //..................
00412   sprintf(tt1,"eq%02dUV",sectID);
00413   sprintf(tt2,"#Sigma from 4-strips, best MIP, plane %02dU+V; eta bin;  #DeltaE MeV",sectID);
00414   h2=new TH2F(tt1,tt2,12,.5,12.5,50,-.1,7.5);
00415   hA[23]=(TH1F*)h2;
00416 
00417   //................
00418   sprintf(tt1,"ca%02d",sectID);
00419   sprintf(tt2,"# UxV candidates per tower, sector=%d; x=spiral=iPhi+60*iEta",sectID);
00420   int mxTw=MaxEtaBins*MaxPhiBins;
00421   hA[24]=new TH1F(tt1,tt2,mxTw,-0.5,mxTw-0.5);    
00422   
00423 
00424   // add histos to the list (if provided)
00425   if(HList) {
00426     for(i=0;i<32;i++) {
00427       if(hA[i]==0) continue;
00428       HList->Add(hA[i]);
00429     }  
00430   }
00431 
00432 }
00433 
00434 
00435 //--------------------------------------------------
00436 //--------------------------------------------------
00437 void EEsmdCal::fillOneTailHisto(char cut, int iEta, int iPhi){
00438   // Fill inclusive spectra for one sector 
00439   int iCut=cut-'a';
00440   assert(iCut>=0 && iCut<kCut);
00441 
00442   int iT=0;
00443   for(iT=0;iT<mxTile;iT++) {
00444     TH1F *h=hT[iCut][iT][iEta][iPhi];
00445     h->Fill(tileAdc[iT][iEta][iPhi]);
00446   }
00447 }
00448 
00449 //--------------------------------------------------
00450 //--------------------------------------------------
00451 void EEsmdCal::fillSmdHisto_a(){
00452   // Fill inclusive spectra for one sector 
00453   int iCut='a'-'a';
00454 
00455   int iuv,istrip;
00456   for(iuv=0;iuv<MaxSmdPlains;iuv++) {
00457     float *adc=smdAdc[iuv];
00458     TH1F **hs=hSs[iCut][iuv]; // single strip spectra
00459     // note, to loop over N-2 strips to get right the sum of pairs 
00460     for(istrip=0;istrip<MaxSmdStrips;istrip++) {
00461       hs[istrip]->Fill(adc[istrip]);
00462     }
00463   }
00464 }
00465 
00466 
00467 //-------------------------------------------------
00468 //-------------------------------------------------
00469 void EEsmdCal::saveHisto(TString fname){
00470   TString outName=fname+".hist.root";
00471   TFile f( outName,"recreate");
00472   assert(f.IsOpen());
00473   printf("%d histos are written  to '%s' ...\n",HList->GetEntries(),outName.Data());
00474   HList->Write();
00475   f.Close();
00476 }