00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00031
00032
00033
00035 #include <Stiostream.h>
00036 #include <math.h>
00037 #include "StEmcHitCollection.h"
00038 #include "St_TableSorter.h"
00039 #include "TArrayF.h"
00040 #include "TArrayL.h"
00041 #include "emc_def.h"
00042 ClassImp(StEmcHitCollection)
00043
00044 StEmcHitCollection::StEmcHitCollection() : St_DataSet("bemc") , StEmcGeom("bemc") {
00045 SetTitle("emc_hit");
00046 mNHit = 0; mEnergySum = 0.0, mEtSum = 0.0;
00047 mModule = 0;
00048 mNumsModule.Set(120); mIndexFirstLast.Set(121);
00049 }
00050
00051 StEmcHitCollection::StEmcHitCollection(const Char_t *Name) : St_DataSet(Name) , StEmcGeom(Name) {
00052
00053 SetTitle("emc_hit");
00054 mNHit = 0; mEnergySum = 0.0, mEtSum = 0.0;
00055 mModule = 0;
00056 mNumsModule.Set(120); mIndexFirstLast.Set(121);
00057 }
00058
00059 StEmcHitCollection::~StEmcHitCollection(){ }
00060
00061 Int_t StEmcHitCollection::ADCtoEnergy(St_emc_hits *emc_hit){
00062
00063 Float_t E;
00064 Int_t i; TString name(GetName());
00065 Int_t n = emc_hit->GetNRows();
00066 emc_hits_st *hit = emc_hit->GetTable();
00067
00068 TString n_ped_h="emc/new/"+name+"/pedestal/ped_"+name+"_h";
00069 TString n_ped ="emc/new/"+name+"/pedestal/ped_"+name;
00070 TString n_slp_h="emc/new/"+name+"/gain/slp_"+name+"_h";
00071 TString n_slp ="emc/new/"+name+"/gain/slp_"+name;
00072
00073 St_DataSetIter local(mEmcCalib);
00074 St_emc_calib_header *m_ped_h = (St_emc_calib_header *) local(n_ped_h);
00075 St_emc_pedestal *m_ped = (St_emc_pedestal *) local(n_ped);
00076 St_emc_calib_header *m_slp_h = (St_emc_calib_header *) local(n_slp_h);
00077 St_emc_adcslope *m_slp = (St_emc_adcslope *) local(n_slp);
00078 if(m_ped_h && m_ped && m_slp_h && m_slp){
00079 emc_calib_header_st *ped_h = m_ped_h->GetTable();
00080 emc_pedestal_st *ped = m_ped->GetTable();
00081 emc_calib_header_st *slp_h = m_slp_h->GetTable();
00082 emc_adcslope_st *slp = m_slp->GetTable();
00083 if(ped_h->det != Detector()){
00084 cout << "***StEmcHitCollection::ADCtoEnergy: Pedestal data is not for this datector: "
00085 << name << endl;
00086 return kStWarn;
00087 }
00088 if(slp_h->det != Detector()){
00089 cout << "***StEmcHitCollection::ADCtoEnergy: Gain data is not for this datector: "
00090 << name << endl;
00091 return kStWarn;
00092 }
00093 if(ped_h->nmodule != NModule() ||
00094 ped_h->neta != NEta() ||
00095 ped_h->nsub != NSub() ||
00096 slp_h->nmodule != NModule() ||
00097 slp_h->neta != NEta() ||
00098 slp_h->nsub != NSub() ){
00099 cout << "***StEmcHitCollection::ADCtoEnergy: Inconsistent Ped or Gain data: "
00100 << name << endl;
00101 return kStWarn;
00102 }
00103
00104 mNHit=0; mEnergySum=0.0; mEtSum =0.0;
00105 switch(slp_h->func){
00106 case 1:
00107
00108 for(i = 0; i<n; i++){
00109 Float_t eta, phi; Int_t id,m,e,s;
00110 m = (Int_t)hit[i].module;
00111 e = (Int_t)hit[i].eta;
00112 s = (Int_t)hit[i].sub;
00113
00114 if(!getId(m, e, s, id)){
00115 getEta(m, e, eta); getPhi(m, s, phi);
00116
00117 if(hit[i].adc>-1){
00118 if(slp[id].p0 > 0.0){
00119 Float_t Et = (hit[i].adc-ped[id].ped)/slp[id].p0;
00120 E = Et * cosh(eta);
00121 if(E<=0.0) E=0.0;
00122 else{
00123 mEnergy[mNHit]=E; mId[mNHit]=id;
00124 mNHit++; mEnergySum+=E; mEtSum+=Et;
00125 }
00126 } else {
00127 cout << "StEmcHitCollection::ADCtoEnergy: Adcslope is less equal 0: "
00128 << name << " : " << id <<endl;
00129 E = -1.0;
00130 }
00131 }else{
00132 E = hit[i].energy;
00133 mEnergy[mNHit]=E; mId[mNHit]=id;
00134 mNHit++; mEnergySum+=E; mEtSum+=E/cosh(eta);
00135 }
00136 } else {
00137 cout<<" Bad index m "<<m<<" e "<<e<<" s "<<s<<endl;
00138 }
00139 }
00140 break;
00141 default:
00142 cout<< "***StEmcHitCollection::ADCtoEnergy: Function has not been implimented yet: "
00143 <<name<<" : "<<slp_h->func<<endl;
00144 return kStWarn;
00145 }
00146 }else{
00147 cout<< "StEmcHitCollection::ADCtoEnergy: Not enough information, use GEANT energy: "
00148 << name << " Nhits " << n << endl;
00149 for(i = 0; i<n; i++){
00150 if(hit[i].energy>0.0){
00151 Float_t eta, phi; Int_t id,m,e,s;
00152 m = (Int_t)hit[i].module;
00153 e = (Int_t)hit[i].eta;
00154 s = (Int_t)hit[i].sub;
00155
00156 if(!getId(m, e, s, id)){
00157 getEta(m, e, eta); getPhi(m, s, phi);
00158 E = hit[i].energy;
00159 if(E>0.0){
00160 mEnergy[mNHit]=E; mId[mNHit]=id;
00161 mNHit++; mEnergySum+=E; mEtSum+=E/cosh(eta);
00162 }
00163 } else {
00164 cout<<" Bad index m "<<m<<" e "<<e<<" s "<<s<<endl;
00165 }
00166 }
00167 }
00168 }
00169 return kStOK;
00170 }
00171
00172 Int_t StEmcHitCollection::fill(St_emc_hits *emc_hits)
00173 {
00174
00175 int i;
00176
00177 Int_t n=emc_hits->GetNRows();
00178 if(n==0) { return kStOK; cout<<" No ADC in "<<GetName()<<endl;}
00179
00180 emc_hits_st *hit = emc_hits->GetTable();
00181 if(hit[0].det<1 || hit[0].det>MAXDET){
00182 cout<< "StEmcHitCollection::Fill: Bad detector# in emc_hits_st : "
00183 << hit[0].det<<" "<< GetName() <<endl;
00184 return kStWarn;
00185 }
00186
00187 mId.Set(n); mEnergy.Set(n);
00188 if(ADCtoEnergy(emc_hits ) != kStOK) return kStWarn;
00189
00190
00191 TArrayF ecopy=mEnergy;
00192 TArrayL idcopy(mNHit);
00193 for(i=0; i< mNHit; i++) {idcopy[i] =(Long_t)mId[i];}
00194
00195 St_TableSorter sort(idcopy.GetArray(), mNHit);
00196 for(i=0; i< mNHit; i++) {
00197 Int_t j=sort.GetIndex(i);
00198 mId[i] = idcopy[j];
00199 mEnergy[i] = ecopy[j];
00200 }
00201
00202
00203 Int_t id, m, e, s, mold;
00204 id = HitId(0);
00205 getBin(id, mold, e, s);
00206 mNumsModule[0] = (Short_t)mold;
00207 mIndexFirstLast[0] = 0;
00208 mIndexFirstLast[1] = mNHit;
00209 mModule = 1;
00210
00211 if(mNHit > 1){
00212 for(i=1; i<mNHit; i++){
00213 id = HitId(i);
00214 getBin(id, m, e, s);
00215 if(m != mold){
00216 mold = m;
00217 mNumsModule[mModule] = (Short_t)mold;
00218 mIndexFirstLast[mModule] = i;
00219 mIndexFirstLast[mModule+1] = mNHit;
00220 ++mModule;
00221 }
00222 }
00223 }
00224
00225
00226
00227 return kStOK;
00228 }
00229
00230 St_emc_hits* StEmcHitCollection::copyToTable(const Char_t *Name){
00231
00232 Int_t i, id, m, e, s, n=0;
00233 emc_hits_st raw;
00234 St_emc_hits *emc_hits = new St_emc_hits((Text_t*)Name, mNHit);
00235 for(i=0; i<mNHit; i++){
00236 if(mEnergy[i]>0.0){
00237 id = (Int_t) mId[i];
00238 getBin(id, m, e, s);
00239 raw.det = Detector();
00240 raw.module = m;
00241 raw.eta = e;
00242 raw.sub = s;
00243 raw.adc = -2;
00244 raw.energy = mEnergy[i];
00245 emc_hits->AddAt(&raw, n);
00246 n++;
00247 }
00248 }
00249 return emc_hits;
00250 }
00251
00252 void StEmcHitCollection::printHits(Int_t n, Int_t start){
00253
00254 Int_t i, m, e, s,id;
00255 cout << endl << GetName() << " : ";
00256 if(mNHit<=0){cout << "No hits" << endl; return;}
00257 else{cout << mNHit << " hits" << " Modules "<<mModule<< endl;}
00258 cout << "Raw ID Module Eta Sub Energy"<<endl;
00259 if(start<0) start=0;
00260 if(n<=0) n=1;
00261 if(start>=mNHit) start=mNHit-1;
00262 if(start+n>=mNHit) n=mNHit-start;
00263
00264 int mold=-999999;
00265 for(i=start; i<start+n; i++){
00266 id = (Int_t) mId[i]; getBin(id, m, e, s);
00267 if(i == start) mold=m;
00268 if(mold != m) {
00269 cout << "-------------------------------"<<endl;
00270 mold = m;
00271 }
00272 cout <<i<<" "<<id<<" "<<m<<" "<<e<<" "<<s<<" "<<mEnergy[i]<<endl;
00273 }
00274
00275 for(i=0; i<mModule; i++){
00276 Int_t jm=mIndexFirstLast[i+1] - mIndexFirstLast[i];
00277 cout<<i<<" #Modules "<<mNumsModule[i]<<" jm "<<jm<<
00278 " First "<<mIndexFirstLast[i]<<" Last "<<mIndexFirstLast[i+1]-1<<endl;
00279 }
00280 }
00281
00282 void StEmcHitCollection::printHitsAll(){
00283
00284 printHits(mNHit, 0);
00285 }
00286
00287 void StEmcHitCollection::Browse(TBrowser *b){
00288
00289 printHits(mNHit, 0);
00290 }
00291
00292 void StEmcHitCollection::printNameTable(){
00293 if(mEmcCalib) printf(" mEmcCalib unzero, name %s \n", mEmcCalib->GetName());
00294 else printf(" <E> Pointer of mEmcCalib is zero ************** \n");
00295 }
