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
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179 #include <assert.h>
00180 #include "StEvent/StTpcHit.h"
00181 #include <algorithm>
00182 #include "StTpcHitMaker.h"
00183
00184 #include "TDataSetIter.h"
00185 #include "StDAQMaker/StDAQReader.h"
00186 #include "TError.h"
00187 #include "string.h"
00188 #include "StEvent.h"
00189 #include "StEvent/StTpcHitCollection.h"
00190 #include "StEvent/StTpcHit.h"
00191 #include "RTS/src/DAQ_TPX/tpxFCF_flags.h"
00192 #include "StTpcRawData.h"
00193 #include "StThreeVectorF.hh"
00194
00195 #include "StDaqLib/TPC/trans_table.hh"
00196 #include "StRtsTable.h"
00197 #include "StDbUtilities/StTpcCoordinateTransform.hh"
00198 #include "StTpcDb/StTpcDb.h"
00199 #include "StDbUtilities/StCoordinates.hh"
00200 #include "StDetectorDbMaker/St_tss_tssparC.h"
00201 #include "StDetectorDbMaker/St_tpcSlewingC.h"
00202 #include "StDetectorDbMaker/St_TpcPadCorrectionC.h"
00203 #include "StDetectorDbMaker/St_tpcPadGainT0C.h"
00204 #include "StDetectorDbMaker/St_tpcAnodeHVavgC.h"
00205 #include "StDetectorDbMaker/St_tpcMaxHitsC.h"
00206 #include "StDetectorDbMaker/StDetectorDbTpcRDOMasks.h"
00207 #include "StDetectorDbMaker/St_tpcPadGainT0C.h"
00208 #include "StDetectorDbMaker/St_tpcPadPlanesC.h"
00209 #include "TFile.h"
00210 #include "TNtuple.h"
00211 #include "TH2.h"
00212 #include "St_tpc_cl.h"
00213 TableClassImpl(St_tpc_cl,tcl_cl);
00214 #include "St_daq_cld.h"
00215 TableClassImpl(St_daq_cld,tcl_cl);
00216 #include "St_daq_sim_cld.h"
00217 TableClassImpl(St_daq_sim_cld,tcl_cl);
00218 #include "St_daq_adc_tb.h"
00219 TableClassImpl(St_daq_adc_tb,daq_adc_tb);
00220 #include "St_daq_sim_adc_tb.h"
00221 TableClassImpl(St_daq_sim_adc_tb,daq_sim_adc_tb);
00222 ClassImp(StTpcHitMaker);
00223 static TNtuple *pulserP = 0;
00224 Float_t StTpcHitMaker::fgDp = .1;
00225 Float_t StTpcHitMaker::fgDt = .2;
00226 Float_t StTpcHitMaker::fgDperp = .1;
00227 static Int_t _debug = 0;
00228
00229
00230
00231 Int_t StTpcHitMaker::Init() {
00232 LOG_INFO << "StTpcHitMaker::Init as\t" << GetName() << endm;
00233 const Char_t *Names[kAll] = {"undef",
00234 "tpc_hits","tpx_hits",
00235 "TpcPulser","TpxPulser",
00236 #if 0
00237 "TpcPadMonitor","TpxPadMonitor",
00238 #endif
00239 "TpcDumpPxls2Nt","TpxDumpPxls2Nt",
00240 "TpxRaw","TpxRaw","TpcAvLaser","TpxAvLaser"};
00241 TString MkName(GetName());
00242 for (Int_t k = 1; k < kAll; k++) {
00243 if (MkName.CompareTo(Names[k],TString::kIgnoreCase) == 0) {kMode = (EMode) k; break;}
00244 }
00245 assert(kMode);
00246 memset(maxHits,0,sizeof(maxHits));
00247 maxBin0Hits = 0;
00248 bin0Hits = 0;
00249 return kStOK ;
00250 }
00251
00252 Int_t StTpcHitMaker::InitRun(Int_t runnumber) {
00253 static Bool_t Done = kFALSE;
00254 SetAttr("minSector",1);
00255 SetAttr("maxSector",24);
00256 SetAttr("minRow",1);
00257 NoRows = St_tpcPadPlanesC::instance()->innerPadRows()+St_tpcPadPlanesC::instance()->outerPadRows();
00258 SetAttr("maxRow",NoRows);
00259 if (! Done) {
00260 if (kMode == kTpxAvLaser || kMode == kTpcAvLaser) {
00261 TFile *f = GetTFile();
00262 if (! f) {gMessMgr->Error() << "with Tpx/Tpc AvLaser you must provide TFile as the 5-th parameter in bfc.C macro" << endm; assert(0);}
00263 f->cd();
00264 enum {NoDim = 3};
00265 const Char_t *NameV[NoDim] = { "row", "pad","time"};
00266 const Int_t nBins[NoDim] = { NoRows, 182, 400};
00267 const Double_t xMin[NoDim] = {0.5 , 0.5, -0.5};
00268 const Double_t xMax[NoDim] = {0.5+NoRows, 182.5, 399.5};
00269 fAvLaser = new THnSparseF *[24];
00270 for (Int_t s = 1; s <= 24; s++) {
00271 fAvLaser[s-1] = new THnSparseF(Form("AvLaser_%02i",s), Form("Averaged laser event for sector %02i",s), NoDim, nBins, xMin, xMax);
00272 fAvLaser[s-1]->CalculateErrors(kTRUE);
00273 for (Int_t i = 0; i < NoDim; i++) {
00274 fAvLaser[s-1]->GetAxis(i)->SetName(NameV[i]);
00275 }
00276 f->Add(fAvLaser[s-1]);
00277 }
00278 }
00279 Done = kTRUE;
00280 }
00281
00282
00283
00284 Int_t maxHitsPerSector = St_tpcMaxHitsC::instance()->maxSectorHits();
00285 Int_t maxBinZeroHits = St_tpcMaxHitsC::instance()->maxBinZeroHits();
00286 Int_t livePads = 0;
00287 Int_t totalPads = 0;
00288 Float_t liveFrac = 1;
00289 for(Int_t sector=1;sector<=24;sector++) {
00290 Int_t liveSecPads = 0;
00291 Int_t totalSecPads = 0;
00292 if (maxHitsPerSector > 0 || maxBinZeroHits > 0) {
00293 for(Int_t row=1;row<=NoRows;row++) {
00294 Int_t numPadsAtRow = St_tpcPadPlanesC::instance()->padsPerRow(row);
00295 totalSecPads += numPadsAtRow;
00296 if (StDetectorDbTpcRDOMasks::instance()->isOn(sector,
00297 StDetectorDbTpcRDOMasks::instance()->rdoForPadrow(row)) &&
00298 St_tpcAnodeHVavgC::instance()->livePadrow(sector,row))
00299 liveSecPads += numPadsAtRow;
00300 }
00301 livePads += liveSecPads;
00302 totalPads += totalSecPads;
00303 }
00304 if (maxHitsPerSector > 0) {
00305 liveFrac = TMath::Max((Float_t) 0.1,
00306 ((Float_t) liveSecPads) / ((Float_t) totalSecPads));
00307 maxHits[sector-1] = (Int_t) (liveFrac * maxHitsPerSector);
00308 if (Debug()) {LOG_INFO << "maxHits in sector " << sector
00309 << " = " << maxHits[sector-1] << endm;}
00310 } else {
00311 maxHits[sector-1] = 0;
00312 if (Debug()) {LOG_INFO << "No maxHits in sector " << sector << endm;}
00313 }
00314 }
00315 if (maxBinZeroHits > 0) {
00316 liveFrac = TMath::Max((Float_t) 0.1,
00317 ((Float_t) livePads) / ((Float_t) totalPads));
00318 maxBin0Hits = (Int_t) (liveFrac * maxBinZeroHits);
00319 if (Debug()) {LOG_INFO << "maxBinZeroHits " << maxBin0Hits << endm;}
00320 } else {
00321 maxBin0Hits = 0;
00322 if (Debug()) {LOG_INFO << "No maxBinZeroHits" << endm;}
00323 }
00324
00325 if (kMode == kTpxAvLaser || kMode == kTpcAvLaser) {
00326 StEvtHddr *header = GetEvtHddr();
00327 if (header) {
00328 TFile *f = GetTFile();
00329 if (! f) {gMessMgr->Error() << "with Tpx/Tpc AvLaser you must provide TFile as the 5-th parameter in bfc.C macro" << endm; assert(0);}
00330 f->cd();
00331 header->Write();
00332 }
00333 }
00334 return kStOK;
00335 }
00336
00337 Int_t StTpcHitMaker::Make() {
00338 Int_t minSector = IAttr("minSector");
00339 Int_t maxSector = IAttr("maxSector");
00340 Int_t minRow = IAttr("minRow");
00341 Int_t maxRow = IAttr("maxRow");
00342
00343 bin0Hits = 0;
00344
00345 for (Int_t sector = minSector; sector <= maxSector; sector++) {
00346 fId = 0;
00347
00348 TString cldadc("cld");
00349 if ( kMode == kTpxRaw || kMode == kTpcRaw ||
00350 kMode == kTpcAvLaser || kMode == kTpxAvLaser) cldadc = "adc";
00351 mQuery = Form("tpx/%s[%i]",cldadc.Data(),sector);
00352 StRtsTable *daqTpcTable = GetNextDaqElement(mQuery);
00353 if (daqTpcTable) {
00354 kReaderType = kStandardTpx;
00355 } else {
00356 mQuery = Form("tpc/legacy[%i]",sector);
00357 daqTpcTable = GetNextDaqElement(mQuery);
00358 if (daqTpcTable) {
00359 kReaderType = kLegacyTpc;
00360 } else {
00361 mQuery = Form("tpx/legacy[%i]",sector);
00362 daqTpcTable = GetNextDaqElement(mQuery);
00363 if (daqTpcTable) {
00364 kReaderType = kLegacyTpx;
00365 }
00366 }
00367 }
00368 Int_t hitsAdded = 0;
00369 while (daqTpcTable) {
00370 assert(Sector() == sector);
00371 Int_t row = NoRows;
00372 fTpc = 0;
00373 if (kReaderType == kLegacyTpx || kReaderType == kLegacyTpc) fTpc = (tpc_t*)*DaqDta()->begin();
00374 else row = daqTpcTable->Row();
00375 if (row >= minRow && row <= maxRow) {
00376 switch (kMode) {
00377 case kTpc:
00378 case kTpx: hitsAdded += UpdateHitCollection(sector); break;
00379 case kTpcPulser:
00380 case kTpxPulser: if (fTpc) DoPulser(sector); break;
00381 #if 0
00382 case kTpcPadMonitor:
00383 case kTpxPadMonitor: if (fTpc) PadMonitor(sector); break;
00384 #endif
00385 case kTpcAvLaser:
00386 case kTpxAvLaser:
00387 if ( fTpc) TpcAvLaser(sector);
00388 else TpxAvLaser(sector);
00389 break;
00390 case kTpcDumpPxls2Nt:
00391 case kTpxDumpPxls2Nt: if (fTpc) DumpPixels2Ntuple(sector); break;
00392 case kTpcRaw:
00393 case kTpxRaw:
00394 if ( fTpc) RawTpcData(sector);
00395 else RawTpxData(sector);
00396 break;
00397 default:
00398 break;
00399 }
00400 }
00401 daqTpcTable = GetNextDaqElement(mQuery);
00402 }
00403 if (maxHits[sector-1] && hitsAdded > maxHits[sector-1]) {
00404 LOG_ERROR << "Too many hits (" << hitsAdded << ") in one sector ("
00405 << sector << "). Skipping event." << endm;
00406 return kStSkip;
00407 }
00408 }
00409 if (maxBin0Hits && bin0Hits > maxBin0Hits) {
00410 LOG_ERROR << "Too many hits (" << bin0Hits
00411 << ") starting at time bin 0. Skipping event." << endm;
00412 return kStSkip;
00413 }
00414 if (kMode == kTpc || kMode == kTpx) {
00415 StEvent *pEvent = dynamic_cast<StEvent *> (GetInputDS("StEvent"));
00416 if (Debug()) {LOG_INFO << "StTpcHitMaker::Make : StEvent has been retrieved " <<pEvent<< endm;}
00417 if (! pEvent) {LOG_INFO << "StTpcHitMaker::Make : StEvent has not been found " << endm; return kStWarn;}
00418 StTpcHitCollection *hitCollection = pEvent->tpcHitCollection();
00419 if (hitCollection) AfterBurner(hitCollection);
00420 }
00421 return kStOK;
00422 }
00423
00424 Int_t StTpcHitMaker::UpdateHitCollection(Int_t sector) {
00425
00426 StEvent *pEvent = dynamic_cast<StEvent *> (GetInputDS("StEvent"));
00427 if (Debug()) {LOG_INFO << "StTpcHitMaker::Make : StEvent has been retrieved " <<pEvent<< endm;}
00428 if (! pEvent) {LOG_INFO << "StTpcHitMaker::Make : StEvent has not been found " << endm; return 0;}
00429 StTpcHitCollection *hitCollection = pEvent->tpcHitCollection();
00430 if ( !hitCollection ) {
00431
00432 hitCollection = new StTpcHitCollection();
00433 pEvent->setTpcHitCollection(hitCollection);
00434 }
00435 Int_t NRows = DaqDta()->GetNRows();
00436 if (NRows <= 0) return 0;
00437 Int_t nhitsBefore = hitCollection->numberOfHits();
00438 Int_t sec = DaqDta()->Sector();
00439 Int_t row = DaqDta()->Row();
00440 if (kReaderType == kLegacyTpc || kReaderType == kLegacyTpx) {
00441 tpc_t *tpc = (tpc_t *) DaqDta()->GetTable();
00442 for (Int_t l = 0; l < NRows; tpc++) {
00443 if ( !tpc->has_clusters ) return 0;
00444 for(Int_t padrow=0;padrow<NoRows;padrow++) {
00445 if (! St_tpcPadGainT0C::instance()->livePadrow(sector,padrow+1)) continue;
00446 tpc_cl *c = &tpc->cl[padrow][0];
00447 Int_t ncounts = tpc->cl_counts[padrow];
00448 for(Int_t j=0;j<ncounts;j++,c++) {
00449 if (! c || ! c->charge) continue;
00450 Int_t iok = hitCollection->addHit(CreateTpcHit(*c,sector,padrow+1));
00451 assert(iok);
00452 }
00453 }
00454 }
00455 } else if (St_tpcPadGainT0C::instance()->livePadrow(sector,row)) {
00456
00457 daq_cld *cld = (daq_cld *) DaqDta()->GetTable();
00458 if (Debug() > 1) {
00459 LOG_INFO << Form("CLD sec %2d: row %2d: clusters: %3d",sec, row, NRows) << endm;
00460 }
00461 for (Int_t l = 0; l < NRows; l++, cld++) {
00462 if (Debug() > 1) {
00463 LOG_INFO << Form(" pad %f[%d:%d], tb %f[%d:%d], cha %d, fla 0x%X",
00464 cld->pad,
00465 cld->p1,
00466 cld->p2,
00467 cld->tb,
00468 cld->t1,
00469 cld->t2,
00470 cld->charge,
00471 cld->flags
00472 ) << endm;
00473 }
00474 if (! cld->pad || ! cld->charge) continue;
00475 if (cld->tb >= __MaxNumberOfTimeBins__) continue;
00476 Int_t iok = hitCollection->addHit(CreateTpcHit(*cld,sector,row));
00477 assert(iok);
00478 }
00479 }
00480 Int_t nhits = hitCollection->numberOfHits() - nhitsBefore;
00481 if (Debug()) {
00482 LOG_INFO << " Total hits in Sector : row " << sector << " : " << row << " = " << nhits << endm;
00483 }
00484 return nhits;
00485 }
00486
00487 StTpcHit *StTpcHitMaker::CreateTpcHit(const tpc_cl &cluster, Int_t sector, Int_t row) {
00488
00489
00490 Float_t pad = cluster.p;
00491 Float_t time = cluster.t;
00492 if (kReaderType == kLegacyTpx) time += 22;
00493 static StTpcCoordinateTransform transform(gStTpcDb);
00494 static StTpcLocalSectorCoordinate local;
00495 static StTpcLocalCoordinate global;
00496 StTpcPadCoordinate padcoord(sector, row, pad, time);
00497 transform(padcoord,local,kFALSE);
00498 transform(local,global);
00499
00500 UInt_t hw = 1;
00501 hw += sector << 4;
00502 hw += row << 9;
00503
00504 Int_t npads = TMath::Abs(cluster.p2 - cluster.p1) + 1;
00505 hw += (npads << 15);
00506
00507 Int_t ntmbk = TMath::Abs(cluster.t2 - cluster.t1) + 1;
00508 hw += (ntmbk << 22);
00509
00510 static StThreeVector<double> hard_coded_errors(fgDp,fgDt,fgDperp);
00511
00512 Double_t gain = (row<=13) ? St_tss_tssparC::instance()->gain_in() : St_tss_tssparC::instance()->gain_out();
00513 Double_t wire_coupling = (row<=13) ? St_tss_tssparC::instance()->wire_coupling_in() : St_tss_tssparC::instance()->wire_coupling_out();
00514 Double_t q = cluster.charge * ((Double_t)St_tss_tssparC::instance()->ave_ion_pot() *
00515 (Double_t)St_tss_tssparC::instance()->scale())/(gain*wire_coupling) ;
00516
00517 StTpcHit *hit = StTpcHitFlag(global.position(),hard_coded_errors,hw,q
00518 , (UChar_t ) 0
00519 , (UShort_t) 0
00520 , (UShort_t) 0
00521 , ++fId
00522 , cluster.p1
00523 , cluster.p2
00524 , cluster.t1
00525 , cluster.t2
00526 , pad
00527 , time
00528 , cluster.charge
00529 , cluster.flags);
00530 if (hit->minTmbk() == 0) bin0Hits++;
00531
00532
00533 return hit;
00534 }
00535
00536 StTpcHit *StTpcHitMaker::CreateTpcHit(const daq_cld &cluster, Int_t sector, Int_t row) {
00537
00538
00539 Float_t pad = cluster.pad;
00540 Float_t time = cluster.tb;
00541
00542 Double_t gain = (row<=13) ? St_tss_tssparC::instance()->gain_in() : St_tss_tssparC::instance()->gain_out();
00543 Double_t wire_coupling = (row<=13) ? St_tss_tssparC::instance()->wire_coupling_in() : St_tss_tssparC::instance()->wire_coupling_out();
00544 Double_t q = cluster.charge * ((Double_t)St_tss_tssparC::instance()->ave_ion_pot() *
00545 (Double_t)St_tss_tssparC::instance()->scale())/(gain*wire_coupling) ;
00546
00547
00548 Double_t freq = gStTpcDb->Electronics()->samplingFrequency();
00549 time = freq * St_tpcSlewingC::instance()->correctedT(row,q,time/freq);
00550
00551 static StTpcCoordinateTransform transform(gStTpcDb);
00552 static StTpcLocalSectorCoordinate local;
00553 static StTpcLocalCoordinate global;
00554 StTpcPadCoordinate padcoord(sector, row, pad, time);
00555 transform(padcoord,local,kFALSE);
00556 transform(local,global);
00557
00558 UInt_t hw = 1;
00559 hw += sector << 4;
00560 hw += row << 9;
00561
00562 Int_t npads = TMath::Abs(cluster.p2 - cluster.p1) + 1;
00563 hw += (npads << 15);
00564
00565 Int_t ntmbk = TMath::Abs(cluster.t2 - cluster.t1) + 1;
00566 hw += (ntmbk << 22);
00567
00568 static StThreeVector<double> hard_coded_errors(fgDp,fgDt,fgDperp);
00569
00570 StTpcHit *hit = StTpcHitFlag(global.position(),hard_coded_errors,hw,q
00571 , (UChar_t ) 0
00572 , (UShort_t) 0
00573 , (UShort_t) 0
00574 , ++fId
00575 , cluster.p1
00576 , cluster.p2
00577 , cluster.t1
00578 , cluster.t2
00579 , pad
00580 , time
00581 , cluster.charge
00582 , cluster.flags);
00583 if (hit->minTmbk() == 0) bin0Hits++;
00584
00585 return hit;
00586 }
00587
00588 void StTpcHitMaker::DoPulser(Int_t sector) {
00589 struct Pulser_t {Float_t sector, row, pad, gain, t0, nnoise, noise, npeak;};
00590 static const Char_t *names = "sector:row:pad:gain:t0:nnoise:noise:npeak";
00591 static Pulser_t Pulser;
00592 if (! pulserP) {
00593 TFile *f = GetTFile();
00594 assert(f);
00595 f->cd();
00596 pulserP = new TNtuple("pulserP","Pulser analysis",names);
00597 }
00598 Int_t r, p, tb, tbmax;
00599 Int_t npeak, nnoise;
00600 if (! fTpc) return;
00601 if (! fTpc->channels_sector) return;
00602 for(Int_t row = 1; row <= NoRows; row++) {
00603 r = row - 1;
00604 if (! fTpc->cl_counts[r]) continue;
00605 for (Int_t pad = 1; pad <= 182; pad++) {
00606 p = pad - 1;
00607 Int_t ncounts = fTpc->counts[r][p];
00608 if (! ncounts) continue;
00609 static UShort_t adc[512];
00610 memset (adc, 0, sizeof(adc));
00611 tbmax = 513;
00612 UShort_t adcmax = 0;
00613 for (Int_t i = 0; i < ncounts; i++) {
00614 tb = fTpc->timebin[r][p][i];
00615 adc[tb] = log8to10_table[fTpc->adc[r][p][i]];
00616 if (adc[tb] > adcmax) {
00617 tbmax = tb;
00618 adcmax = adc[tb];
00619 }
00620 }
00621 if (tbmax < 2 || tbmax > 504) continue;
00622 npeak = nnoise = 0;
00623 Int_t i1s = TMath::Max( 0, tbmax - 2);
00624 Int_t i2s = TMath::Min(511, tbmax + 7);
00625 Int_t i1 = TMath::Max(0 ,i1s - 20);
00626 Int_t i2 = TMath::Min(511,i2s + 20);
00627 Double_t peak = 0;
00628 Double_t noise = 0;
00629 Double_t t0 = 0;
00630 for (Int_t i = i1; i <= i2; i++) {
00631 if (i >= i1s && i <= i2s) continue;
00632 nnoise++;
00633 noise += adc[i];
00634 }
00635 if (nnoise) noise /= nnoise;
00636 for (Int_t i = i1s; i <= i2s; i++) {
00637 npeak++;
00638 peak += adc[i] - noise;
00639 t0 += i*(adc[i] - noise);
00640 }
00641 if (peak <= 0) continue;
00642 t0 /= peak;
00643 Pulser.sector = sector;
00644 Pulser.row = row;
00645 Pulser.pad = pad;
00646 Pulser.gain = peak;
00647 Pulser.t0 = t0;
00648 Pulser.nnoise = nnoise;
00649 Pulser.noise = noise;
00650 Pulser.npeak = npeak;
00651 pulserP->Fill(&Pulser.sector);
00652 }
00653 }
00654 }
00655 #if 0
00656
00657 void StTpcHitMaker::PadMonitor(Int_t sector) {
00658 static TH2F *padMon[24][NoRows];
00659 static TH2F *pcl[24][NoRows];
00660 static Bool_t first = kTRUE;
00661 if (! fTpc) return;
00662 if (first) {
00663 first = kFALSE;
00664 for (Int_t s = 1; s <= 24; s++)
00665 for (Int_t r = 1; r <= NoRows; r++) {
00666 padMon[s-1][r-1] = new TH2F(Form("padS%02iR%02i",s,r),Form("Pad monitor for sector = %i and row = %i",s,r),
00667 512,0,512,182,1,183);
00668 pcl[s-1][r-1] = new TH2F(Form("clS%02iR%02i",s,r),Form("Cluster monitor for sector = %i and row = %i",s,r),
00669 512,0,512,182,1,183);
00670 }
00671 } else {
00672 for (Int_t s = 1; s <= 24; s++)
00673 for (Int_t r = 1; r <= NoRows; r++) {
00674 padMon[s-1][r-1]->Reset();
00675 pcl[s-1][r-1]->Reset();
00676 }
00677 }
00678 Int_t nhits = 0;
00679 Int_t npixels = 0;
00680
00681 Int_t s = sector - 1;
00682 if (Debug()) PrintSpecial(sector);
00683 if ( fTpc->has_clusters ) {
00684 for(Int_t r=0;r<NoRows;r++) {
00685 tpc_cl *c = &fTpc->cl[r][0];
00686 Int_t ncounts = fTpc->cl_counts[r];
00687 for(Int_t j=0;j<ncounts;j++,c++) {
00688 Float_t pad = c->p;
00689 Float_t time = c->t;
00690 Float_t q = c->charge;
00691 pcl[s][r]->Fill(time,pad,q);
00692 nhits++;
00693 }
00694 }
00695 LOG_INFO << " Total hits in Sector : " << sector << " = " << nhits << endm;
00696 }
00697
00698 for(Int_t r = 0; r < NoRows; r++) {
00699 for (Int_t pad = 1; pad <= 182; pad++) {
00700 Int_t p = pad - 1;
00701 Int_t ncounts = fTpc->counts[r][p];
00702 if (! ncounts) continue;
00703 for (Int_t i = 0; i < ncounts; i++) {
00704 Int_t tb = fTpc->timebin[r][p][i];
00705 Float_t adc = log8to10_table[fTpc->adc[r][p][i]];
00706 padMon[s][r]->Fill(tb,pad,adc);
00707 npixels++;
00708 }
00709 }
00710 }
00711 LOG_INFO << " Total pixels in Sector : " << sector << " = " << npixels << endm;
00712 }
00713 #endif
00714
00715 void StTpcHitMaker::TpcAvLaser(Int_t sector) {
00716 if (! fTpc) return;
00717 Int_t npixels = 0;
00718 struct pixl_t {
00719 Double_t sector, row, pad, time;
00720 };
00721 if (fAvLaser[sector-1]->GetNbins() > 500000) {
00722 THnSparseF *hnew = CompressTHn(fAvLaser[sector-1]);
00723 GetTFile()->Remove(fAvLaser[sector-1]);
00724 delete fAvLaser[sector-1];
00725 fAvLaser[sector-1] = hnew;
00726 GetTFile()->Add(fAvLaser[sector-1]);
00727 }
00728 pixl_t pixel;
00729 pixel.sector = sector;
00730 for(Int_t r = 0; r < NoRows; r++) {
00731 pixel.row = r+1;
00732 for (Int_t pad = 1; pad <= 182; pad++) {
00733 pixel.pad = pad;
00734 Int_t p = pad - 1;
00735 Double_t gain = St_tpcPadGainT0C::instance()->Gain(pixel.sector,pixel.row,pixel.pad);
00736 if (gain <= 0) continue;
00737 Int_t ncounts = fTpc->counts[r][p];
00738 if (! ncounts) continue;
00739 for (Int_t i = 0; i < ncounts; i++) {
00740 pixel.time = fTpc->timebin[r][p][i];
00741 Double_t adc = log8to10_table[fTpc->adc[r][p][i]];
00742 fAvLaser[sector-1]->Fill(&pixel.row,gain*adc);
00743 npixels++;
00744 }
00745 }
00746 }
00747 LOG_INFO << " Total pixels in Sector : " << sector << " = " << npixels << endm;
00748 }
00749
00750 void StTpcHitMaker::TpxAvLaser(Int_t sector) {
00751 assert(fAvLaser[sector-1]);
00752 if (fAvLaser[sector-1]->GetNbins() > 1000000) {
00753 THnSparseF *hnew = CompressTHn(fAvLaser[sector-1]);
00754 GetTFile()->Remove(fAvLaser[sector-1]);
00755 delete fAvLaser[sector-1];
00756 fAvLaser[sector-1] = hnew;
00757 GetTFile()->Add(fAvLaser[sector-1]);
00758 }
00759 Int_t r=Row() ;
00760 if(r==0) return;
00761 r-- ;
00762 Int_t p = Pad() - 1 ;
00763 if (p < 0 || p >= St_tpcPadPlanesC::instance()->padsPerRow(r+1)) return;
00764 struct pixl_t {
00765 Double_t sector, row, pad, time;
00766 };
00767 pixl_t pixel;
00768 pixel.sector = sector;
00769 pixel.row = r+1;
00770 pixel.pad = p+1;
00771 TGenericTable::iterator iword = DaqDta()->begin();
00772 for (;iword != DaqDta()->end();++iword) {
00773 daq_adc_tb &daqadc = (*(daq_adc_tb *)*iword);
00774 pixel.time = daqadc.tb;
00775 Double_t adc = daqadc.adc;
00776 if (adc < 6) continue;
00777 fAvLaser[sector-1]->Fill(&pixel.row,adc);
00778 }
00779 }
00780
00781 void StTpcHitMaker::DumpPixels2Ntuple(Int_t sector) {
00782 struct BPoint_t {
00783 Float_t sector, row, pad, tb, adc, ped, t0, peak;
00784 };
00785 static const Char_t *BName = "sector:row:pad:tb:adc:ped:t0:peak";
00786 static TNtuple *adcP = 0;
00787 if (! adcP) {
00788 assert(GetTFile());
00789 GetTFile()->cd();
00790 adcP = new TNtuple("adcP","Pulser ADC",BName);
00791 }
00792 static BPoint_t P;
00793 if (! fTpc) return;
00794 Int_t r, p, tb, tbmax;
00795
00796 for(Int_t row = 1; row <= NoRows; row++) {
00797 r = row - 1;
00798 for (Int_t pad = 1; pad <= 182; pad++) {
00799 p = pad - 1;
00800 Int_t ncounts = fTpc->counts[r][p];
00801 if (! ncounts) continue;
00802 static UShort_t adc[512];
00803 memset (adc, 0, sizeof(adc));
00804 tbmax = 513;
00805 UShort_t adcmax = 0;
00806 for (Int_t i = 0; i < ncounts; i++) {
00807 tb = fTpc->timebin[r][p][i];
00808 adc[tb] = log8to10_table[fTpc->adc[r][p][i]];
00809 if (adc[tb] > adcmax) {
00810 tbmax = tb;
00811 adcmax = adc[tb];
00812 }
00813 }
00814 if (tbmax < 2 || tbmax > 504) continue;
00815 Int_t npeak = 0, nped = 0;
00816 Int_t i1s = TMath::Max( 0, tbmax - 2);
00817 Int_t i2s = TMath::Min(511, tbmax + 7);
00818 Int_t i1 = TMath::Max(0 ,i1s - 20);
00819 Int_t i2 = TMath::Min(511,i2s + 20);
00820 Double_t peak = 0;
00821 Double_t ped = 0;
00822 Double_t t0 = 0;
00823 for (Int_t i = i1; i <= i2; i++) {
00824 if (i >= i1s && i <= i2s) continue;
00825 nped++;
00826 ped += adc[i];
00827 }
00828 if (nped) ped /= nped;
00829 for (Int_t i = i1s; i <= i2s; i++) {
00830 npeak++;
00831 peak += adc[i] - ped;
00832 t0 += i*(adc[i] - ped);
00833 }
00834 if (peak <= 0) continue;
00835 t0 /= peak;
00836 i1 = (Int_t) TMath::Max(0.,t0 - 20);
00837 i2 = (Int_t) TMath::Min(511., t0 + 80);
00838 for (Int_t i = i1; i <= i2; i++) {
00839 P.sector = sector;
00840 P.row = row;
00841 P.pad = pad;
00842 P.tb = i - t0;
00843 P.adc = adc[i];
00844 P.ped = ped;
00845 P.t0 = t0;
00846 P.peak = peak;
00847 adcP->Fill(&P.sector);
00848 }
00849 }
00850 }
00851 }
00852
00853 void StTpcHitMaker::PrintSpecial(Int_t sector) {
00854
00855
00856 Int_t r,p,t ;
00857 UInt_t adc = 0;
00858 UChar_t val ;
00859 if (! fTpc) return;
00860 if(fTpc->mode==0) {
00861 UInt_t tot_pix = 0 ;
00862 UInt_t cl_count = 0 ;
00863 Int_t i ;
00864
00865 for(r=0;r<NoRows;r++) {
00866 for(p=0;p<182;p++) {
00867 for(t=0;t<fTpc->counts[r][p];t++) {
00868 val = fTpc->adc[r][p][t] ;
00869 Int_t vali = log8to10_table[val];
00870 adc += val ;
00871 if(val) tot_pix++ ;
00872 if (Debug() > 1) {
00873 Int_t timebin = fTpc->timebin[r][p][t] ;
00874 printf("%d %d %d %d %d\n",sector,r+1,p+1,timebin,vali) ;
00875 }
00876 }
00877 }
00878
00879 if(fTpc->has_clusters) {
00880 cl_count += fTpc->cl_counts[r] ;
00881 }
00882 if (Debug() > 1) {
00883 if(fTpc->has_clusters) {
00884 for(i=0;i<fTpc->cl_counts[r];i++) {
00885 tpc_cl *c = &fTpc->cl[r][i] ;
00886
00887 printf("%d %d %f %f %d %d %d %d %d %d\n",
00888 sector,r+1,c->p,c->t,c->charge,c->flags,c->p1,c->p2,c->t1,c->t2) ;
00889 }
00890 }
00891 }
00892 }
00893 LOG_INFO << Form("TPC: Sector %d: occupancy %3d %%, charge %d, pixels %u, clusters %d",sector,
00894 (int)(100.0 *((double)fTpc->channels_sector/(double)fTpc->max_channels_sector)),
00895 adc,tot_pix,cl_count) << endm;
00896 }
00897 }
00898
00899 StTpcDigitalSector *StTpcHitMaker::GetDigitalSector(Int_t sector) {
00900 TDataSet *event = GetData("Event");
00901 StTpcRawData *data = 0;
00902 if (! event) {
00903 data = new StTpcRawData(24);
00904 event = new TObjectSet("Event", data);
00905 AddData(event);
00906 } else data = (StTpcRawData *) event->GetObject();
00907 assert(data);
00908 StTpcDigitalSector *digitalSector = data->GetSector(sector);
00909 if (! digitalSector) {
00910 digitalSector = new StTpcDigitalSector();
00911 data->setSector(sector,digitalSector);
00912 }
00913 return digitalSector;
00914 }
00915
00916 Int_t StTpcHitMaker::RawTpxData(Int_t sector) {
00917 memset(ADCs, 0, sizeof(ADCs));
00918 memset(IDTs, 0, sizeof(IDTs));
00919 StTpcDigitalSector *digitalSector = 0;
00920 Int_t some_data = 0;
00921 Int_t r_old = -1;
00922 Int_t p_old = -1;
00923 Int_t Total_data = 0;
00924 Int_t r=Row() ;
00925 if(r==0) return 0 ;
00926 r-- ;
00927 Int_t p = Pad() - 1 ;
00928 if (p < 0 || p >= St_tpcPadPlanesC::instance()->padsPerRow(r+1)) return 0;
00929 if (r_old != r || p_old != p) {
00930 if (some_data) {
00931 Total_data += some_data;
00932 some_data = 0;
00933 if (! digitalSector) digitalSector = GetDigitalSector(sector);
00934 Int_t ntbold = digitalSector->numberOfTimeBins(r_old+1,p_old+1);
00935 if (ntbold) {
00936 LOG_INFO << "digitalSector " << sector
00937 << " already has " << ntbold << " at row/pad " << r_old+1 << "/" << p_old+1 << endm;
00938 }
00939 digitalSector->putTimeAdc(r_old+1,p_old+1,ADCs,IDTs);
00940 memset(ADCs, 0, sizeof(ADCs));
00941 memset(IDTs, 0, sizeof(IDTs));
00942 }
00943 r_old = r;
00944 p_old = p;
00945 }
00946 TGenericTable::iterator iword = DaqDta()->begin();
00947 for (;iword != DaqDta()->end();++iword) {
00948 daq_adc_tb &daqadc = (*(daq_adc_tb *)*iword);
00949 Int_t tb = daqadc.tb;
00950 Int_t adc = daqadc.adc;
00951 ADCs[tb] = adc;
00952 IDTs[tb] = 65535;
00953 some_data++ ;
00954 }
00955
00956 if (some_data) {
00957 Total_data += some_data;
00958 some_data = 0;
00959 if (! digitalSector) digitalSector = GetDigitalSector(sector);
00960 Int_t ntbold = digitalSector->numberOfTimeBins(r_old+1,p_old+1);
00961 if (ntbold) {
00962 LOG_INFO << "digitalSector " << sector
00963 << " already has " << ntbold << " at row/pad " << r_old+1 << "/" << p_old+1 << endm;
00964 }
00965 digitalSector->putTimeAdc(r_old+1,p_old+1,ADCs,IDTs);
00966 memset(ADCs, 0, sizeof(ADCs));
00967 memset(IDTs, 0, sizeof(IDTs));
00968 }
00969 return Total_data;
00970 }
00971
00972 Int_t StTpcHitMaker::RawTpcData(Int_t sector) {
00973 if (! fTpc) return 0;
00974 memset(ADCs, 0, sizeof(ADCs));
00975 memset(IDTs, 0, sizeof(IDTs));
00976 StTpcDigitalSector *digitalSector = 0;
00977 Int_t Total_data = 0;
00978 for (Int_t row = 1; row <= NoRows; row++) {
00979 Int_t r = row - 1;
00980 if (! digitalSector) digitalSector = GetDigitalSector(sector);
00981 for (Int_t pad = 1; pad <= digitalSector->numberOfPadsAtRow(row); pad++) {
00982 Int_t p = pad - 1;
00983 memset(ADCs, 0, sizeof(ADCs));
00984 memset(IDTs, 0, sizeof(IDTs));
00985 Int_t ncounts = fTpc->counts[r][p];
00986 if (! ncounts) continue;
00987 for (Int_t i = 0; i < ncounts; i++) {
00988 Int_t tb = fTpc->timebin[r][p][i];
00989 ADCs[tb] = log8to10_table[fTpc->adc[r][p][i]];
00990 IDTs[tb] = 65535;
00991 Total_data++;
00992 }
00993 Int_t ntbold = digitalSector->numberOfTimeBins(row,pad);
00994 if (ntbold) {
00995 LOG_INFO << "digitalSector " << sector
00996 << " already has " << ntbold << " at row/pad " << row << "/" << pad << endm;
00997 }
00998 digitalSector->putTimeAdc(row,pad,ADCs,IDTs);
00999 }
01000 }
01001 if (Total_data) {
01002 LOG_INFO << "Read " << Total_data << " pixels from Sector " << sector << endm;
01003 }
01004 return Total_data;
01005 }
01006
01007 #ifdef __MAKE_NTUPLE__
01008 static TNtuple *tup = 0;
01009 struct TpcHitPair_t {
01010 Float_t sec, row,
01011 qK, padK, tbkK, padKmn, padKmx, tbkKmn, tbkKmx, IdTK, QAK,
01012 qL, padL, tbkL, padLmn, padLmx, tbkLmn, tbkLmx, IdTL, QAL,
01013 padOv, tbkOv;
01014 };
01015 static const Char_t *vTpcHitMRPair = "sec:row:"
01016 "qK:padK:tbkK:padKmn:padKmx:tbkKmn:tbkKmx:IdTK:QAK:"
01017 "qL:padL:tbkL:padLmn:padLmx:tbkLmn:tbkLmx:IdTL:QAL:"
01018 "padOv:tbkOv";
01019 #endif
01020
01021 Bool_t TpcHitLess(const StTpcHit *lhs, const StTpcHit *rhs) {
01022 return (200*lhs->timeBucket() + lhs->pad()) < (200*rhs->timeBucket() + rhs->pad());
01023 };
01024
01025 void StTpcHitMaker::AfterBurner(StTpcHitCollection *TpcHitCollection) {
01026 static Float_t padDiff = 2.5, timeBucketDiff = 5.0;
01027 static StTpcCoordinateTransform transform(gStTpcDb);
01028 static StTpcLocalSectorCoordinate local;
01029 static StTpcLocalCoordinate global;
01030 if (! TpcHitCollection) return;
01031 #ifdef __MAKE_NTUPLE__
01032 if (! tup) {
01033 if (StChain::GetChain()->GetTFile()) {
01034 StChain::GetChain()->GetTFile()->cd();
01035 tup = new TNtuple("HitT","Cluster Pair Info",vTpcHitMRPair);
01036 }
01037 }
01038 TpcHitPair_t pairC;
01039 #endif
01040 UInt_t numberOfSectors = TpcHitCollection->numberOfSectors();
01041 for (UInt_t sec = 1; sec <= numberOfSectors; sec++) {
01042 StTpcSectorHitCollection* sectorCollection = TpcHitCollection->sector(sec-1);
01043 if (sectorCollection) {
01044 UInt_t numberOfPadrows = sectorCollection->numberOfPadrows();
01045 for (UInt_t row = 1; row <= numberOfPadrows; row++) {
01046 StTpcPadrowHitCollection *rowCollection = TpcHitCollection->sector(sec-1)->padrow(row-1);
01047 if (rowCollection) {
01048 UInt_t NoHits = rowCollection->hits().size();
01049 if (NoHits < 2) continue;
01050 sort(rowCollection->hits().begin(),
01051 rowCollection->hits().end(),
01052 TpcHitLess);
01053
01054 Int_t merged = 0;
01055 for (UInt_t k = 0; k < NoHits; k++) {
01056 StTpcHit* kHit = TpcHitCollection->sector(sec-1)->padrow(row-1)->hits().at(k);
01057 if (_debug) {cout << "k " << k; kHit->Print();}
01058 if (kHit->flag()) continue;
01059 #ifdef __MAKE_NTUPLE__
01060 pairC.sec = sec;
01061 pairC.row = row;
01062 pairC.qK = kHit->charge();
01063 pairC.padK = kHit->pad();
01064 pairC.tbkK = kHit->timeBucket();
01065 pairC.padKmn = kHit->minPad();
01066 pairC.padKmx = kHit->maxPad();
01067 pairC.tbkKmn = kHit->minTmbk();
01068 pairC.tbkKmx = kHit->maxTmbk();
01069 pairC.IdTK = kHit->idTruth();
01070 pairC.QAK = kHit->qaTruth();
01071 #endif
01072 for (UInt_t l = 0; l < NoHits; l++) {
01073 if (k == l) continue;
01074 StTpcHit* lHit = TpcHitCollection->sector(sec-1)->padrow(row-1)->hits().at(l);
01075 if (_debug) {cout << "l " << l; lHit->Print();}
01076 if (lHit->flag()) continue;
01077
01078 Int_t padOverlap = TMath::Min(kHit->maxPad(),lHit->maxPad())
01079 - TMath::Max(kHit->minPad(),lHit->minPad());
01080 if (padOverlap < 0) continue;
01081 Int_t tmbkOverlap = TMath::Min(kHit->maxTmbk(),lHit->maxTmbk())
01082 - TMath::Max(kHit->minTmbk(),lHit->minTmbk());
01083 if (tmbkOverlap < 0) continue;
01084 #ifdef __MAKE_NTUPLE__
01085 if (tup) {
01086 pairC.qL = lHit->charge();
01087 pairC.padL = lHit->pad();
01088 pairC.tbkL = lHit->timeBucket();
01089 pairC.padLmn = lHit->minPad();
01090 pairC.padLmx = lHit->maxPad();
01091 pairC.tbkLmn = lHit->minTmbk();
01092 pairC.tbkLmx = lHit->maxTmbk();
01093 pairC.IdTL = lHit->idTruth();
01094 pairC.QAL = lHit->qaTruth();
01095 pairC.padOv = padOverlap;
01096 pairC.tbkOv = tmbkOverlap;
01097 tup->Fill(&pairC.sec);
01098 }
01099 #endif
01100
01101 if (TMath::Abs(kHit->pad() - lHit->pad()) > padDiff ||
01102 TMath::Abs(kHit->timeBucket() - lHit->timeBucket()) > timeBucketDiff) continue;
01103 #ifdef FCF_CHOPPED
01104 UShort_t flag = lHit->flag() | FCF_CHOPPED;
01105 #else
01106 UShort_t flag = lHit->flag() | 0x080;
01107 #endif
01108 lHit->setFlag(flag);
01109 if (_debug) {
01110 cout << "mk" << k; kHit->Print();
01111 cout << "ml" << l; lHit->Print();
01112 }
01113 Float_t q = lHit->charge() + kHit->charge();
01114 Float_t adc = lHit->adc() + kHit->adc();
01115 Float_t pad = (lHit->charge()*lHit->pad() + kHit->charge()*kHit->pad() )/q;
01116 Float_t timeBucket = (lHit->charge()*lHit->timeBucket() + kHit->charge()*kHit->timeBucket())/q;
01117 Short_t minPad = TMath::Min(lHit->minPad(), kHit->minPad());
01118 Short_t maxPad = TMath::Max(lHit->maxPad(), kHit->maxPad());
01119 Short_t minTmbk = TMath::Min(lHit->minTmbk(), kHit->minTmbk());
01120 Short_t maxTmbk = TMath::Max(lHit->maxTmbk(), kHit->maxTmbk());
01121 Int_t IdT = lHit->idTruth();
01122 Double_t QA = lHit->charge()*lHit->qaTruth();
01123 if (IdT == kHit->idTruth()) {
01124 QA += kHit->charge()*kHit->qaTruth();
01125 } else {
01126 if (lHit->charge()*lHit->qaTruth() < kHit->charge()*kHit->qaTruth()) {
01127 QA = kHit->charge()*kHit->qaTruth();
01128 IdT = kHit->idTruth();
01129 }
01130 }
01131 QA = QA/q;
01132 kHit->setIdTruth(IdT, TMath::Nint(QA));
01133 kHit->setCharge(q);
01134 kHit->setExtends(pad,timeBucket,minPad,maxPad,minTmbk,maxTmbk);
01135 kHit->setAdc(adc);
01136 StTpcPadCoordinate padcoord(sec, row, pad, timeBucket);
01137 transform(padcoord,local,kFALSE);
01138 transform(local,global);
01139 kHit->setPosition(global.position());
01140 if (_debug) {
01141 cout << "m " << k; kHit->Print();
01142 }
01143 merged++;
01144 }
01145 }
01146 #ifdef __CORRECT_S_SHAPE__
01147
01148 for (UInt_t k = 0; k < NoHits; k++) {
01149 StTpcHit* kHit = TpcHitCollection->sector(sec-1)->padrow(row-1)->hits().at(k);
01150 if (kHit->flag()) continue;
01151 Double_t pad = kHit->pad();
01152 Double_t timeBucket = kHit->timeBucket();
01153 Int_t io = 1;
01154 if (row > 13) io = 2;
01155 Int_t np = kHit->padsInHit();
01156 pad += St_TpcPadCorrectionC::instance()->GetCorrection(pad,io,np,0);
01157 StTpcPadCoordinate padcoord(sec, row, pad, timeBucket);
01158 transform(padcoord,local,kFALSE);
01159 transform(local,global);
01160 kHit->setPosition(global.position());
01161 }
01162 #endif
01163 }
01164 }
01165 }
01166 }
01167 return;
01168 };
01169
01170 StTpcHit* StTpcHitMaker::StTpcHitFlag(const StThreeVectorF& p,
01171 const StThreeVectorF& e,
01172 UInt_t hw, float q, UChar_t c,
01173 UShort_t idTruth, UShort_t quality,
01174 UShort_t id,
01175 Short_t mnpad, Short_t mxpad, Short_t mntmbk,
01176 Short_t mxtmbk, Float_t cl_x, Float_t cl_t, UShort_t adc,
01177 UShort_t flag) {
01178
01179 StTpcHit* hit = new StTpcHit(p,e,hw,q,c,idTruth,quality,id,mnpad,mxpad,mntmbk,mxtmbk,cl_x,cl_t,adc);
01180
01181
01182 if ( mntmbk<0 || mxtmbk<0 || mntmbk>32000 || mxtmbk>32000
01183 || mnpad <0 || mxpad <0 || mnpad >32000 || mxpad >32000
01184 || mxpad-mnpad > 100
01185 || (Float_t) mntmbk>cl_t || (Float_t) mxtmbk<cl_t
01186 || (Float_t) mnpad >cl_x || (Float_t) mxpad <cl_x
01187 ) flag |= FCF_SANITY;
01188 hit->setFlag(flag);
01189 return hit;
01190 }
01191
01192 THnSparseF *StTpcHitMaker::CompressTHn(THnSparseF *hist, Double_t compress) {
01193 if (! hist) return 0;
01194 Int_t nd = hist->GetNdimensions();
01195 Int_t *nbins = new Int_t[nd];
01196 for (Int_t i = 0; i < nd; i++) nbins[i] = hist->GetAxis(i)->GetNbins();
01197 THnSparseF *hnew = new THnSparseF(hist->GetName(),hist->GetTitle(),nd, nbins, 0, 0);
01198 hnew->CalculateErrors(kTRUE);
01199 delete [] nbins;
01200 for (Int_t i = 0; i < nd; i++) {
01201 TAxis *ax = hist->GetAxis(i);
01202 if (ax->IsVariableBinSize()) hnew->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
01203 else hnew->GetAxis(i)->Set(ax->GetNbins(), ax->GetXmin(), ax->GetXmax());
01204 }
01205 Int_t *bins = new Int_t[nd];
01206 Double_t *x = new Double_t[nd];
01207 Long64_t N = hist->GetNbins(); cout << hist->GetName() << " has " << N << " bins before compression." << endl;
01208 Double_t max = -1;
01209 for (Long64_t i = 0; i < N; ++i) {
01210 Double_t cont = hist->GetBinContent(i, bins);
01211 if (cont > max) max = cont;
01212 }
01213 for (Long64_t i = 0; i < N; ++i) {
01214 Double_t cont = hist->GetBinContent(i, bins);
01215 if (cont < max/compress) continue;
01216
01217 for (Int_t d = 0; d < nd; ++d) {x[d] = hist->GetAxis(d)->GetBinCenter(bins[d]);}
01218 hnew->Fill(x,cont);
01219 }
01220 delete [] bins;
01221 delete [] x;
01222 cout << hnew->GetName() << " has " << hnew->GetNbins() << " bins after compression." << endl;
01223 return hnew;
01224 }
