// Class for Evaluation and Validation of the ALTRO Tail Cancelation Filter //
// (TCF) parameters out of TPC Raw data //
// //
-// Author: Stefan Rossegger //
+// Author: Stefan Rossegger, Simon Feigl //
// //
///////////////////////////////////////////////////////////////////////////////
#include <TStyle.h>
#include <TMinuit.h>
#include <TH1F.h>
+#include <TH2F.h>
+#include <AliSysInfo.h>
#include <TMath.h>
#include <TNtuple.h>
#include <TEntryList.h>
-
#include "AliRawReaderRoot.h"
+#include "AliRawHLTManager.h"
#include "AliTPCRawStream.h"
+#include "AliTPCRawStreamV3.h"
#include "AliTPCROC.h"
#include "AliTPCAltroEmulator.h"
+#include "AliTPCmapper.h"
+#include <fstream>
+
ClassImp(AliTPCCalibTCF)
AliTPCCalibTCF::AliTPCCalibTCF() :
TNamed(),
- fGateWidth(100),
+ fGateWidth(50),
fSample(900),
- fPulseLength(500),
+ fPulseLength(400),
fLowPulseLim(30),
- fUpPulseLim(1000),
- fRMSLim(4.)
+ fUpPulseLim(900),
+ fRMSLim(1.0),
+ fRatioIntLim(2)
+
{
//
// AliTPCCalibTCF standard constructor
}
//_____________________________________________________________________________
-AliTPCCalibTCF::AliTPCCalibTCF(Int_t gateWidth, Int_t sample, Int_t pulseLength, Int_t lowPulseLim, Int_t upPulseLim, Double_t rmsLim) :
+AliTPCCalibTCF::AliTPCCalibTCF(Int_t gateWidth, Int_t sample, Int_t pulseLength, Int_t lowPulseLim, Int_t upPulseLim, Double_t rmsLim, Double_t ratioIntLim) :
TNamed(),
fGateWidth(gateWidth),
fSample(sample),
fPulseLength(pulseLength),
fLowPulseLim(lowPulseLim),
fUpPulseLim(upPulseLim),
- fRMSLim(rmsLim)
+ fRMSLim(rmsLim),
+ fRatioIntLim(ratioIntLim)
{
//
// AliTPCCalibTCF constructor with specific (non-standard) thresholds
fPulseLength(tcf.fPulseLength),
fLowPulseLim(tcf.fLowPulseLim),
fUpPulseLim(tcf.fUpPulseLim),
- fRMSLim(tcf.fRMSLim)
+ fRMSLim(tcf.fRMSLim),
+ fRatioIntLim(tcf.fRatioIntLim)
{
//
// AliTPCCalibTCF copy constructor
//
}
+
+//_____________________________________________________________________________
+void AliTPCCalibTCF::ProcessRawFileV3(const char *nameRawFile, const char *nameFileOut) {
+ //
+ // New RCU data format!: Standard middle of 2009
+ //
+ // Loops over all events within one RawData file and collects proper pulses
+ // (according to given tresholds) per pad
+ // Histograms per pad are stored in 'nameFileOut'
+ //
+
+ AliRawReader *rawReader = AliRawReader::Create(nameRawFile);
+ if (!rawReader) {
+ printf("Could not create a raw reader for %s\n",nameRawFile);
+ return;
+ }
+
+ rawReader->RewindEvents(); // just to make sure
+
+ rawReader->Select("TPC");
+
+ if (!rawReader->NextEvent()) {
+ printf("no events found in %s\n",nameRawFile);
+ return;
+ }
+
+ // TPC stream reader
+ AliTPCRawStreamV3 rawStream(rawReader);
+
+ Int_t ievent=0;
+ do {
+ AliSysInfo::AddStamp(Form("start_event_%d",ievent), ievent,-1,-1);
+ printf("Reading next event ... Nr: %d\n",ievent);
+ // Start the basic data extraction
+ ProcessRawEventV3(rawReader, &rawStream, nameFileOut);
+ AliSysInfo::AddStamp(Form("end_event_%d",ievent), ievent,-1,-1);
+ ievent++;
+
+ } while (rawReader->NextEvent());
+
+ rawReader->~AliRawReader();
+
+}
+
+
//_____________________________________________________________________________
-void AliTPCCalibTCF::ProcessRawFile(const char *nameRawFile, const char *nameFileOut) {
+void AliTPCCalibTCF::ProcessRawFile(const char *nameRawFile, const char *nameFileOut, bool bUseHLTOUT) {
//
// Loops over all events within one RawData file and collects proper pulses
// (according to given tresholds) per pad
// Histograms per pad are stored in 'nameFileOut'
//
+ // create the data reader
AliRawReader *rawReader = new AliRawReaderRoot(nameRawFile);
- rawReader->Reset();
+ if (!rawReader) {
+ return;
+ }
+
+ // create HLT reader for redirection of TPC data from HLTOUT to TPC reconstruction
+ AliRawReader *hltReader=AliRawHLTManager::AliRawHLTManager::CreateRawReaderHLT(rawReader, "TPC");
- while ( rawReader->NextEvent() ){ // loop
- printf("Reading next event ...");
- AliTPCRawStream rawStream(rawReader);
- rawReader->Select("TPC");
- ProcessRawEvent(&rawStream, nameFileOut);
+ // now choose the data source
+ if (bUseHLTOUT) rawReader=hltReader;
+
+ // rawReader->Reset();
+ rawReader->RewindEvents();
+
+ if (!rawReader->NextEvent()) {
+ printf("no events found in %s\n",nameRawFile);
+ return;
}
+ Int_t ievent=0;
+ do {
+ AliSysInfo::AddStamp(Form("start_event_%d",ievent), ievent,-1,-1);
+ printf("Reading next event ... Nr: %d\n",ievent);
+ AliTPCRawStream *rawStream = new AliTPCRawStream(rawReader);
+ rawReader->Select("TPC");
+ ProcessRawEvent(rawStream, nameFileOut);
+ delete rawStream;
+ AliSysInfo::AddStamp(Form("end_event_%d",ievent), ievent,-1,-1);
+ ievent++;
+ } while (rawReader->NextEvent());
+
rawReader->~AliRawReader();
}
+//_____________________________________________________________________________
+void AliTPCCalibTCF::ProcessRawEventV3( AliRawReader *rawReader, AliTPCRawStreamV3 *rawStream, const char *nameFileOut) {
+ //
+ // New RCU data format!: Standard middle of 2009
+ //
+ // Extracts proper pulses (according the given tresholds) within one event
+ // and accumulates them into one histogram per pad. All histograms are
+ // saved in the file 'nameFileOut'.
+ // The first bins of the histograms contain the following information:
+ // bin 1: Number of accumulated pulses
+ // bin 2;3;4: Sector; Row; Pad;
+ //
+
+ TFile fileOut(nameFileOut,"UPDATE");
+ fileOut.cd();
+
+ TH1I *tempHis = new TH1I("tempHis","tempHis",fSample,fGateWidth,fSample+fGateWidth);
+ TH1I *tempRMSHis = new TH1I("tempRMSHis","tempRMSHis",2000,0,2000);
+
+ // loop over the data in this event
+
+ while (rawStream->NextDDL() ) {
+
+ Int_t ddl = rawReader->GetDDLID();
+
+ while (rawStream->NextChannel() ) {
+
+ while (rawStream->NextBunch() ) {
+
+ Int_t t0 = rawStream->GetStartTimeBin();
+ Int_t bl = rawStream->GetBunchLength();
+
+ if (bl<fSample+fGateWidth) continue;
+
+ Int_t sector = rawStream->GetSector();
+ Int_t row = rawStream->GetRow();
+ Int_t pad = rawStream->GetPad();
+
+ UShort_t *signals=(UShort_t*)rawStream->GetSignals();
+ if (!signals) continue;
+
+ // Write to temporary histogramm
+ for (Int_t i=0;i<bl;++i) {
+ UShort_t time=t0-i;
+ UShort_t signal=signals[i];
+ if ( (fGateWidth<time) && (time<=fSample+fGateWidth) ) {
+ tempHis->SetBinContent(time-fGateWidth,signal);
+ }
+ }
+
+ // calculation of the pulse properties and comparison to thresholds settings
+
+ Int_t max = (Int_t)tempHis->GetMaximum(FLT_MAX);
+ Int_t maxpos = tempHis->GetMaximumBin();
+
+ Int_t first = (Int_t)TMath::Max(maxpos-10, 0);
+ Int_t last = TMath::Min((Int_t)maxpos+fPulseLength-10, fSample+fGateWidth);
+
+ // simple baseline substraction ? better one needed ? (pedestalsubstr.?)
+ // and RMS calculation with timebins before the pulse and at the end of
+ // the signal
+ for (Int_t ipos = 0; ipos<6; ipos++) {
+ // before the pulse
+ tempRMSHis->Fill(tempHis->GetBinContent(first+ipos));
+ }
+ for (Int_t ipos = 0; ipos<20; ipos++) {
+ // at the end to get rid of pulses with serious baseline fluctuations
+ tempRMSHis->Fill(tempHis->GetBinContent(last-ipos));
+ }
+
+ Double_t baseline = tempRMSHis->GetMean();
+ Double_t rms = tempRMSHis->GetRMS();
+ tempRMSHis->Reset();
+
+ Double_t lowLim = fLowPulseLim+baseline;
+ Double_t upLim = fUpPulseLim+baseline;
+
+ // get rid of pulses which contain gate signal and/or too much noise
+ // with the help of ratio of integrals
+ Double_t intHist = 0;
+ Double_t intPulse = 0;
+ Double_t binValue;
+ for(Int_t ipos=first; ipos<=last; ipos++) {
+ binValue = TMath::Abs(tempHis->GetBinContent(ipos) - baseline);
+ intHist += binValue;
+ if(ipos>=first+5 && ipos<=first+15) {intPulse += binValue;}
+ }
+
+ // gets rid of high frequency noise:
+ // calculating ratio (value one to the right of maximum)/(maximum)
+ // has to be >= 0.1; if maximum==0 set ratio to 0.1
+ Double_t maxCorr = max - baseline;
+ Double_t binRatio = 0.1;
+ if(TMath::Abs(maxCorr)>1e-5) {
+ binRatio = (tempHis->GetBinContent(maxpos+1) - baseline) / maxCorr;
+ }
+
+ // Decision if found pulse is a proper one according to given tresholds
+ if (max>lowLim && max<upLim && !((last-first)<fPulseLength) && rms<fRMSLim && (intHist/intPulse)<fRatioIntLim &&intPulse>10&& (binRatio >= 0.1) ) {
+
+ // 1D histogramm for mean pulse per pad
+ char hname[100];
+ snprintf(hname,100,"sec%drow%dpad%d",sector,row,pad);
+
+ TH1F *his = (TH1F*)fileOut.Get(hname);
+
+ if (!his ) { // new entry (pulse in new pad found)
+
+ his = new TH1F(hname,hname, fPulseLength+5, 0, fPulseLength+5);
+ his->SetBinContent(1,1); // pulse counter (1st pulse)
+ his->SetBinContent(2,sector); // sector
+ his->SetBinContent(3,row); // row
+ his->SetBinContent(4,pad); // pad
+
+ for (Int_t ipos=0; ipos<last-first; ipos++){
+ Int_t signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline);
+ his->SetBinContent(ipos+5,signal);
+ }
+ his->Write(hname);
+ printf("new %s: Signal %d at bin %d \n", hname, max-(Int_t)baseline, maxpos+fGateWidth);
+
+ } else { // adding pulse to existing histogram (pad already found)
+
+ his->AddBinContent(1,1); // pulse counter for each pad
+ for (Int_t ipos=0; ipos<last-first; ipos++){
+ Int_t signal= (Int_t)(tempHis->GetBinContent(ipos+first)-baseline);
+ his->AddBinContent(ipos+5,signal);
+ }
+ printf("adding ... %s: Signal %d at bin %d \n", hname, max-(Int_t)baseline, maxpos+fGateWidth);
+ his->Write(hname,kOverwrite);
+ }
+
+
+ // 2D histogramm for pulse spread within a DDL (normalized to one)
+ char hname2d[100];
+ snprintf(hname2d,100,"2Dhisto_ddl%d",ddl);
+ TH2F *his2d = (TH2F*)fileOut.Get(hname2d);
+ if (!his2d ) { // new entry (ddl was not seen before)
+
+ his2d = new TH2F(hname2d,hname2d, fPulseLength, 0., (Double_t)fPulseLength, 50,-0.02,0.02);
+ for (Int_t ipos=0; ipos<last-first; ipos++){
+ Double_t signal = tempHis->GetBinContent(ipos+first)-baseline;
+ if (TMath::Abs(signal/maxCorr)>1e-10) // zero bins are biased
+ his2d->Fill(ipos,signal/maxCorr);
+ }
+ his2d->Write(hname2d);
+ printf("new %s: \n", hname2d);
+ } else { // adding pulse to existing histogram
+
+ for (Int_t ipos=0; ipos<last-first; ipos++){
+ Double_t signal= tempHis->GetBinContent(ipos+first)-baseline;
+ if (TMath::Abs(signal/maxCorr)>1e-10) // zero bins are biased
+ his2d->Fill(ipos,signal/maxCorr);
+ }
+ his2d->Write(hname2d,kOverwrite);
+ }
+
+ tempHis->Reset();
+
+ } // pulse stored
+
+ } // bunch loop
+ }// channel loop
+ } // ddl loop
+
+ tempHis->~TH1I();
+ tempRMSHis->~TH1I();
+ printf("Finished to read event ... \n");
+ fileOut.Close();
+
+}
+
+
//_____________________________________________________________________________
void AliTPCCalibTCF::ProcessRawEvent(AliTPCRawStream *rawStream, const char *nameFileOut) {
//
// bin 2;3;4: Sector; Row; Pad;
//
+ rawStream->Reset();
+
Int_t sector = rawStream->GetSector();
Int_t row = rawStream->GetRow();
+
+ Int_t prevSec = 999999;
+ Int_t prevRow = 999999;
+ Int_t prevPad = 999999;
Int_t prevTime = 999999;
- Int_t prevPad = 999999;
TFile fileOut(nameFileOut,"UPDATE");
fileOut.cd();
TH1I *tempHis = new TH1I("tempHis","tempHis",fSample+fGateWidth,fGateWidth,fSample+fGateWidth);
TH1I *tempRMSHis = new TH1I("tempRMSHis","tempRMSHis",2000,0,2000);
- rawStream->SetOldRCUFormat(1);
+ // printf("raw next: %d\n",rawStream->Next());
while (rawStream->Next()) {
if (rawStream->IsNewRow()){
sector = rawStream->GetSector();
row = rawStream->GetRow();
+ // if (sector!=prevSec) AliSysInfo::AddStamp(Form("sector_%d_row_%d",sector,row), -1,sector,row);
}
Int_t pad = rawStream->GetPad();
Int_t time = rawStream->GetTime();
Int_t signal = rawStream->GetSignal();
- if (!rawStream->IsNewPad()) { // Reading signal from one Pad
+ //printf("%d: t:%d sig:%d\n",pad,time,signal);
+
+ // Reading signal from one Pad
+ if (!rawStream->IsNewPad()) {
+
+ // this pad always gave a useless signal, probably induced by the supply
+ // voltage of the gate signal (date:2008-Aug-07)
+ if(sector==51 && row==95 && pad==0) {
+ continue;
+ }
+
+ // only process pulses of pads with correct address
+ if(sector<0 || sector+1 > Int_t(AliTPCROC::Instance()->GetNSector())) {
+ continue;
+ }
+ if(row<0 || row+1 > Int_t(AliTPCROC::Instance()->GetNRows(sector))) {
+ continue;
+ }
+ if(pad<0 || pad+1 > Int_t(AliTPCROC::Instance()->GetNPads(sector,row))) {
+ continue;
+ }
+
if (time>prevTime) {
- printf("Wrong time: %d %d\n",rawStream->GetTime(),prevTime);
- rawStream->Dump();
+ //printf("Wrong time: %d %d\n",rawStream->GetTime(),prevTime);
+ continue;
} else {
// still the same pad, save signal to temporary histogram
- if (time<=fSample+fGateWidth && time>fGateWidth) {
+ if ( (time<=fSample+fGateWidth) && (time>=fGateWidth)) {
tempHis->SetBinContent(time,signal);
}
- }
+ }
+
} else {
+
// complete pulse found and stored into tempHis, now calculation
- // of it's properties and comparison to given thresholds
-
+ // of the properties and comparison to given thresholds
+
Int_t max = (Int_t)tempHis->GetMaximum(FLT_MAX);
Int_t maxpos = tempHis->GetMaximumBin();
Int_t first = (Int_t)TMath::Max(maxpos-10, 0);
- Int_t last = TMath::Min((Int_t)maxpos+fPulseLength-10, fSample);
+ Int_t last = TMath::Min((Int_t)maxpos+fPulseLength-10, fSample+fGateWidth);
// simple baseline substraction ? better one needed ? (pedestalsubstr.?)
// and RMS calculation with timebins before the pulse and at the end of
for (Int_t ipos = 0; ipos<6; ipos++) {
// before the pulse
tempRMSHis->Fill(tempHis->GetBinContent(first+ipos));
+ }
+ for (Int_t ipos = 0; ipos<20; ipos++) {
// at the end to get rid of pulses with serious baseline fluctuations
tempRMSHis->Fill(tempHis->GetBinContent(last-ipos));
}
+
Double_t baseline = tempRMSHis->GetMean();
Double_t rms = tempRMSHis->GetRMS();
tempRMSHis->Reset();
Double_t lowLim = fLowPulseLim+baseline;
Double_t upLim = fUpPulseLim+baseline;
+ // get rid of pulses which contain gate signal and/or too much noise
+ // with the help of ratio of integrals
+ Double_t intHist = 0;
+ Double_t intPulse = 0;
+ Double_t binValue;
+ for(Int_t ipos=first; ipos<=last; ipos++) {
+ binValue = TMath::Abs(tempHis->GetBinContent(ipos) - baseline);
+ intHist += binValue;
+ if(ipos>=first+5 && ipos<=first+25) {intPulse += binValue;}
+ }
+
+ // gets rid of high frequency noise:
+ // calculating ratio (value one to the right of maximum)/(maximum)
+ // has to be >= 0.1; if maximum==0 set ratio to 0.1
+ Double_t maxCorr = max - baseline;
+ Double_t binRatio = 0.1;
+ if(TMath::Abs(maxCorr)>1e-5) {
+ binRatio = (tempHis->GetBinContent(maxpos+1) - baseline) / maxCorr;
+ }
+
// Decision if found pulse is a proper one according to given tresholds
- if (max>lowLim && max<upLim && !((last-first)<fPulseLength) && rms<fRMSLim){
+ if (max>lowLim && max<upLim && !((last-first)<fPulseLength) && rms<fRMSLim && (intHist/intPulse)<fRatioIntLim && (binRatio >= 0.1) ) {
char hname[100];
- sprintf(hname,"sec%drow%dpad%d",sector,row,prevPad);
+ snprintf(hname,100,"sec%drow%dpad%d",prevSec,prevRow,prevPad);
TH1F *his = (TH1F*)fileOut.Get(hname);
if (!his ) { // new entry (pulse in new pad found)
his = new TH1F(hname,hname, fPulseLength+4, 0, fPulseLength+4);
- his->SetBinContent(1,1); // pulse counter (1st pulse)
- his->SetBinContent(2,sector); // sector
- his->SetBinContent(3,row); // row
- his->SetBinContent(4,prevPad); // pad
+ his->SetBinContent(1,1); // pulse counter (1st pulse)
+ his->SetBinContent(2,prevSec); // sector
+ his->SetBinContent(3,prevRow); // row
+ his->SetBinContent(4,prevPad); // pad
for (Int_t ipos=0; ipos<last-first; ipos++){
- Int_t signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline);
+ signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline);
his->SetBinContent(ipos+5,signal);
}
his->Write(hname);
his->AddBinContent(1,1); // pulse counter for each pad
for (Int_t ipos=0; ipos<last-first; ipos++){
- Int_t signal= (Int_t)(tempHis->GetBinContent(ipos+first)-baseline);
+ signal= (Int_t)(tempHis->GetBinContent(ipos+first)-baseline);
his->AddBinContent(ipos+5,signal);
}
printf("adding ... %s: Signal %d at bin %d \n", hname, max-(Int_t)baseline, maxpos+fGateWidth);
tempHis->Reset();
}
prevTime = time;
+ prevSec = sector;
+ prevRow = row;
prevPad = pad;
}
TIter next( fileIn.GetListOfKeys() );
char nameFileOut[100];
- sprintf(nameFileOut,"Sec-%s",nameFileIn);
+ snprintf(nameFileOut,100,"Sec-%s",nameFileIn);
TFile fileOut(nameFileOut,"RECREATE");
fileOut.cd();
while ( (key=(TKey*)next()) ) {
iHist++;
+ TString name(key->GetName());
+ if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl
+
+ hisPad = (TH1F*)fileIn.Get(name.Data()); // copy object to memory
- hisPad = (TH1F*)fileIn.Get(key->GetName()); // copy object to memory
Int_t pulseLength = hisPad->GetNbinsX() -4;
// -4 because first four timebins contain pad specific informations
Int_t npulse = (Int_t)hisPad->GetBinContent(1);
Int_t sector = (Int_t)hisPad->GetBinContent(2);
char hname[100];
- sprintf(hname,"sector%d",sector);
+ snprintf(hname,100,"sector%d",sector);
TH1F *his = (TH1F*)fileOut.Get(hname);
if (!his ) { // new histogram (new sector)
for (Int_t ipos=0; ipos<pulseLength; ipos++){
his->AddBinContent(ipos+5,hisPad->GetBinContent(ipos+5));
}
- his->Write(hname,kOverwrite);
+ his->Write(hname,kOverwrite);
}
if (iHist%500==0) {
printf("merging status: \t %d pads out of %d \n",iHist, nHist);
}
}
+
printf("merging done ...\n");
fileIn.Close();
fileOut.Close();
- // calculate TCF parameters on averaged pulse per Sector
- AnalyzeRootFile(nameFileOut);
-
}
//____________________________________________________________________________
-void AliTPCCalibTCF::AnalyzeRootFile(const char *nameFileIn, Int_t minNumPulse) {
+void AliTPCCalibTCF::AnalyzeRootFile(const char *nameFileIn, Int_t minNumPulse, Int_t histStart, Int_t histEnd) {
//
// This function takes a prepeared root file (accumulated histograms: output
// of process function) and performs an analysis (fit and equalization) in
TIter next( fileIn.GetListOfKeys() );
char nameFileOut[100];
- sprintf(nameFileOut,"TCFparam-%s",nameFileIn);
+ snprintf(nameFileOut,100,"TCF-%s",nameFileIn);
TFile fileOut(nameFileOut,"RECREATE");
fileOut.cd();
Int_t nHist = fileIn.GetNkeys();
Int_t iHist = 0; // counter for print of analysis-status
- while (key = (TKey *) next()) { // loop over histograms
-
- printf("Analyze histogramm %d out of %d\n",++iHist,nHist);
- hisIn = (TH1F*)fileIn.Get(key->GetName()); // copy object to memory
+ while ((key = (TKey *) next())) { // loop over histograms
+ ++iHist;
+ if(iHist < histStart || iHist > histEnd) {continue;}
+ TString name(key->GetName());
+ if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl
+
+ hisIn = (TH1F*)fileIn.Get(key->GetName()); // copy object to memory
+
Int_t numPulse = (Int_t)hisIn->GetBinContent(1);
if ( numPulse >= minNumPulse ) {
-
+ printf("Analyze histogram %d out of %d\n",iHist,nHist);
Double_t* coefP = new Double_t[3];
Double_t* coefZ = new Double_t[3];
for(Int_t i = 0; i < 3; i++){
}
coefP->~Double_t();
coefZ->~Double_t();
+ } else {
+ printf("Skip histogram %d out of %d | not enough accumulated pulses\n",iHist,nHist);
}
-
+
}
fileIn.Close();
//____________________________________________________________________________
-Int_t AliTPCCalibTCF::AnalyzePulse(TH1F *hisIn, Double_t *coefZ, Double_t *coefP) {
+Int_t AliTPCCalibTCF::AnalyzePulse(TH1F * const hisIn, Double_t *coefZ, Double_t *coefP) {
//
// Performs the analysis on one specific pulse (histogram) by means of fitting
// the pulse and equalization of the pulseheight. The found TCF parameters
//
Int_t pulseLength = hisIn->GetNbinsX() -4;
- // -1 because the first four timebins usually contain pad specific informations
+ // -4 because the first four timebins usually contain pad specific informations
Int_t npulse = (Int_t)hisIn->GetBinContent(1);
Int_t sector = (Int_t)hisIn->GetBinContent(2);
Int_t row = (Int_t)hisIn->GetBinContent(3);
if (fitOk) {
// calculates the 3rd set (remaining 2 PZ values) in order to restore the
// original height of the pulse
- Equalization(dataTuple, coefZ, coefP);
-
+ Int_t equOk = Equalization(dataTuple, coefZ, coefP);
+ if (!equOk) {
+ Error("FindFit", "Pulse equalisation procedure failed - pulse abandoned ");
+ printf("in Sector %d | Row %d | Pad %d |", sector, row, pad);
+ printf(" Npulses: %d \n\n", npulse);
+ coefP[2] = 0; coefZ[2] = 0;
+ dataTuple->~TNtuple();
+ return 0;
+ }
printf("Calculated TCF parameters for: \n");
printf("Sector %d | Row %d | Pad %d |", sector, row, pad);
printf(" Npulses: %d \n", npulse);
//____________________________________________________________________________
-void AliTPCCalibTCF::TestTCFonRootFile(const char *nameFileIn, const char *nameFileTCF, Int_t plotFlag, Int_t lowKey, Int_t upKey)
+void AliTPCCalibTCF::TestTCFonRootFile(const char *nameFileIn, const char *nameFileTCF, Int_t minNumPulse, Int_t plotFlag, Int_t lowKey, Int_t upKey)
{
//
// Performs quality parameters evaluation of the calculated TCF parameters in
}
char nameFileOut[100];
- sprintf(nameFileOut,"Quality_%s_AT_%s",nameFileTCF, nameFileIn);
+ snprintf(nameFileOut,100,"Quality_%s_AT_%s",nameFileTCF, nameFileIn);
TFile fileOut(nameFileOut,"RECREATE");
TNtuple *qualityTuple = new TNtuple("TCFquality","TCF quality Values","sec:row:pad:npulse:heightDev:areaRed:widthRed:undershot:maxUndershot");
Int_t iHist = 0;
for(Int_t i=0;i<lowKey-1;i++){++iHist; key = (TKey *) next();}
- while (key = (TKey *) next()) { // loop over saved histograms
+ while ((key = (TKey *) next())) { // loop over saved histograms
// loading pulse to memory;
+ TString name(key->GetName());
+ if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl
+
printf("validating pulse %d out of %d\n",++iHist,nHist);
hisIn = (TH1F*)fileIn.Get(key->GetName());
+
// find the correct TCF parameter according to the his infos (first 4 bins)
Int_t nPulse = FindCorTCFparam(hisIn, nameFileTCF, coefZ, coefP);
- if (nPulse) { // doing the TCF quality analysis
+ if (nPulse>=minNumPulse) { // doing the TCF quality analysis
Double_t *quVal = GetQualityOfTCF(hisIn,coefZ,coefP, plotFlag);
Int_t sector = (Int_t)hisIn->GetBinContent(2);
Int_t row = (Int_t)hisIn->GetBinContent(3);
//_____________________________________________________________________________
-void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameFileOut, const char *nameFileTCF, Int_t plotFlag) {
+void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameFileOut, const char *nameFileTCF, Int_t minNumPulse, Int_t plotFlag, bool bUseHLTOUT) {
//
// Performs quality parameters evaluation of the calculated TCF parameters in
// the file 'nameFileTCF' for every proper pulse (according to given thresholds)
// and test the found TCF parameters on them ...
//
+
+ // create the data reader
AliRawReader *rawReader = new AliRawReaderRoot(nameRawFile);
- rawReader->Reset();
-
+ if (!rawReader) {
+ return;
+ }
+
+ // create HLT reader for redirection of TPC data from HLTOUT to TPC reconstruction
+ AliRawReader *hltReader=AliRawHLTManager::AliRawHLTManager::CreateRawReaderHLT(rawReader, "TPC");
+
+ // now choose the data source
+ if (bUseHLTOUT) rawReader=hltReader;
+
+ // rawReader->Reset();
+ rawReader->RewindEvents();
+
+ if (!rawReader->NextEvent()) {
+ printf("no events found in %s\n",nameRawFile);
+ return;
+ }
+
Double_t* coefP = new Double_t[3];
Double_t* coefZ = new Double_t[3];
for(Int_t i = 0; i < 3; i++){
coefZ[i] = 0;
}
- while ( rawReader->NextEvent() ){
+ Int_t ievent = 0;
+
+ TH1I *tempHis = new TH1I("tempHis","tempHis",fSample+fGateWidth,fGateWidth,fSample+fGateWidth);
+ TH1I *tempRMSHis = new TH1I("tempRMSHis","tempRMSHis",2000,0,2000);
+
+ TFile fileOut(nameFileOut,"UPDATE"); // Quality Parameters storage
+ TNtuple *qualityTuple = (TNtuple*)fileOut.Get("TCFquality");
+ if (!qualityTuple) { // no entry in file
+ qualityTuple = new TNtuple("TCFquality","TCF quality Values","sec:row:pad:npulse:heightDev:areaRed:widthRed:undershot:maxUndershot:pulseRMS");
+ }
+
+ do {
- printf("Reading next event...");
- AliTPCRawStream rawStream(rawReader);
+ printf("Reading next event ... Nr:%d\n",ievent);
+ AliTPCRawStream *rawStream = new AliTPCRawStream(rawReader);
rawReader->Select("TPC");
+ ievent++;
- Int_t sector = rawStream.GetSector();
- Int_t row = rawStream.GetRow();
+ Int_t sector = rawStream->GetSector();
+ Int_t row = rawStream->GetRow();
+
+ Int_t prevSec = 999999;
+ Int_t prevRow = 999999;
+ Int_t prevPad = 999999;
Int_t prevTime = 999999;
- Int_t prevPad = 999999;
-
- TH1I *tempHis = new TH1I("tempHis","tempHis",fSample+fGateWidth,fGateWidth,fSample+fGateWidth);
- TH1I *tempRMSHis = new TH1I("tempRMSHis","tempRMSHis",2000,0,2000);
-
- rawStream.SetOldRCUFormat(1);
-
- TFile fileOut(nameFileOut,"UPDATE"); // Quality Parameters storage
- TNtuple *qualityTuple = (TNtuple*)fileOut.Get("TCFquality");
- if (!qualityTuple) { // no entry in file
- qualityTuple = new TNtuple("TCFquality","TCF quality Values","sec:row:pad:npulse:heightDev:areaRed:widthRed:undershot:maxUndershot:pulseRMS");
- }
- while (rawStream.Next()) {
+ while (rawStream->Next()) {
- if (rawStream.IsNewRow()){
- sector = rawStream.GetSector();
- row = rawStream.GetRow();
+ if (rawStream->IsNewRow()){
+ sector = rawStream->GetSector();
+ row = rawStream->GetRow();
}
- Int_t pad = rawStream.GetPad();
- Int_t time = rawStream.GetTime();
- Int_t signal = rawStream.GetSignal();
+ Int_t pad = rawStream->GetPad();
+ Int_t time = rawStream->GetTime();
+ Int_t signal = rawStream->GetSignal();
- if (!rawStream.IsNewPad()) { // Reading signal from one Pad
+ if (!rawStream->IsNewPad()) { // Reading signal from one Pad
+
+ // this pad always gave a useless signal, probably induced by the supply
+ // voltage of the gate signal (date:2008-Aug-07)
+ if(sector==51 && row==95 && pad==0) {
+ continue;
+ }
+
+ // only process pulses of pads with correct address
+ if(sector<0 || sector+1 > Int_t(AliTPCROC::Instance()->GetNSector())) {
+ continue;
+ }
+ if(row<0 || row+1 > Int_t(AliTPCROC::Instance()->GetNRows(sector))) {
+ continue;
+ }
+ if(pad<0 || pad+1 > Int_t(AliTPCROC::Instance()->GetNPads(sector,row))) {
+ continue;
+ }
+
if (time>prevTime) {
- printf("Wrong time: %d %d\n",rawStream.GetTime(),prevTime);
- rawStream.Dump();
+ // printf("Wrong time: %d %d\n",rawStream->GetTime(),prevTime);
+ continue;
} else {
+ // still the same pad, save signal to temporary histogram
if (time<=fSample+fGateWidth && time>fGateWidth) {
tempHis->SetBinContent(time,signal);
}
Int_t maxpos = tempHis->GetMaximumBin();
Int_t first = (Int_t)TMath::Max(maxpos-10, 0);
- Int_t last = TMath::Min((Int_t)maxpos+fPulseLength-10, fSample);
+ Int_t last = TMath::Min((Int_t)maxpos+fPulseLength-10, fSample+fGateWidth);
// simple baseline substraction ? better one needed ? (pedestalsubstr.?)
for (Int_t ipos = 0; ipos<6; ipos++) {
// before the pulse
tempRMSHis->Fill(tempHis->GetBinContent(first+ipos));
+ }
+ for (Int_t ipos = 0; ipos<20; ipos++) {
// at the end to get rid of pulses with serious baseline fluctuations
- tempRMSHis->Fill(tempHis->GetBinContent(last-ipos));
+ tempRMSHis->Fill(tempHis->GetBinContent(last-ipos));
}
Double_t baseline = tempRMSHis->GetMean();
Double_t rms = tempRMSHis->GetRMS();
Double_t lowLim = fLowPulseLim+baseline;
Double_t upLim = fUpPulseLim+baseline;
+ // get rid of pulses which contain gate signal and/or too much noise
+ // with the help of ratio of integrals
+ Double_t intHist = 0;
+ Double_t intPulse = 0;
+ Double_t binValue;
+ for(Int_t ipos=first; ipos<=last; ipos++) {
+ binValue = TMath::Abs(tempHis->GetBinContent(ipos) - baseline);
+ intHist += binValue;
+ if(ipos>=first+5 && ipos<=first+15) {intPulse += binValue;}
+ }
+
+ // gets rid of high frequency noise:
+ // calculating ratio (value one to the right of maximum)/(maximum)
+ // has to be >= 0.1; if maximum==0 set ratio to 0.1
+ Double_t maxCorr = max - baseline;
+ Double_t binRatio = 0.1;
+ if(TMath::Abs(maxCorr) > 1e-5 ) {
+ binRatio = (tempHis->GetBinContent(maxpos+1) - baseline) / maxCorr;
+ }
+
+
// Decision if found pulse is a proper one according to given tresholds
- if (max>lowLim && max<upLim && !((last-first)<fPulseLength) && rms<fRMSLim){
+ if (max>lowLim && max<upLim && !((last-first)<fPulseLength) && rms<fRMSLim && intHist/intPulse<fRatioIntLim && (binRatio >= 0.1) ){
// note:
// assuming that lowLim is higher than the pedestal value!
char hname[100];
- sprintf(hname,"sec%drow%dpad%d",sector,row,prevPad);
+ snprintf(hname,100,"sec%drow%dpad%d",prevSec,prevRow,prevPad);
TH1F *his = new TH1F(hname,hname, fPulseLength+4, 0, fPulseLength+4);
his->SetBinContent(1,1); // pulse counter (1st pulse)
- his->SetBinContent(2,sector); // sector
- his->SetBinContent(3,row); // row
- his->SetBinContent(4,prevPad); // pad
+ his->SetBinContent(2,prevSec); // sector
+ his->SetBinContent(3,prevRow); // row
+ his->SetBinContent(4,prevPad); // pad
+
for (Int_t ipos=0; ipos<last-first; ipos++){
- Int_t signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline);
+ signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline);
his->SetBinContent(ipos+5,signal);
}
// (first 4 bins)
Int_t nPulse = FindCorTCFparam(his, nameFileTCF, coefZ, coefP);
- if (nPulse) { // Parameters found - doing the TCF quality analysis
+ if (nPulse>=minNumPulse) { // Parameters found - doing the TCF quality analysis
Double_t *quVal = GetQualityOfTCF(his,coefZ,coefP, plotFlag);
qualityTuple->Fill(sector,row,pad,nPulse,quVal[0],quVal[1],quVal[2],quVal[3],quVal[4],quVal[5]);
quVal->~Double_t();
tempHis->Reset();
}
prevTime = time;
+ prevSec = sector;
+ prevRow = row;
prevPad = pad;
+
}
- tempHis->~TH1I();
- tempRMSHis->~TH1I();
-
- printf("Finished to read event - close output file ... \n");
-
- fileOut.cd();
- qualityTuple->Write("TCFquality",kOverwrite);
- fileOut.Close();
+ printf("Finished to read event ... \n");
+ delete rawStream;
- } // event loop
+ } while (rawReader->NextEvent()); // event loop
+ printf("Finished to read file - close output file ... \n");
+
+ fileOut.cd();
+ qualityTuple->Write("TCFquality",kOverwrite);
+ fileOut.Close();
+
+ tempHis->~TH1I();
+ tempRMSHis->~TH1I();
coefP->~Double_t();
coefZ->~Double_t();
}
+//____________________________________________________________________________
+TH2F *AliTPCCalibTCF::PlotOccupSummary2Dhist(const char *nameFileIn, Int_t side) {
+ //
+ // Plots the number of summed pulses per pad on a given TPC side
+ // 'nameFileIn': root-file created with the Process function
+ //
+
+ TFile fileIn(nameFileIn,"READ");
+ TH1F *his;
+ TKey *key;
+ TIter next(fileIn.GetListOfKeys());
+
+ TH2F * his2D = new TH2F("his2D","his2D", 250,-250,250,250,-250,250);
+
+ AliTPCROC * roc = AliTPCROC::Instance();
+
+ Int_t nHist=fileIn.GetNkeys();
+ if (!nHist) { return 0; }
+
+ Int_t iHist = 0;
+ Float_t xyz[3];
+
+ Int_t binx = 0;
+ Int_t biny = 0;
+
+ Int_t npulse = 0;
+ Int_t sec = 0;
+ Int_t row = 0;
+ Int_t pad = 0;
+
+ while ((key = (TKey *) next())) { // loop over histograms within the file
+ iHist++;
+
+ TString name(key->GetName());
+ if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl
+
+ his = (TH1F*)fileIn.Get(key->GetName()); // copy object to memory
+
+ npulse = (Int_t)his->GetBinContent(1);
+ sec = (Int_t)his->GetBinContent(2);
+ row = (Int_t)his->GetBinContent(3);
+ pad = (Int_t)his->GetBinContent(4);
+
+ if ( (side==0) && (sec%36>=18) ) continue;
+ if ( (side>0) && (sec%36<18) ) continue;
+
+ if ( (row<0) && (pad<0) ) { // row and pad are equal to -1, then -> summed pulses per sector
+ // fill all pad with this values
+ for (UInt_t rowi=0; rowi<roc->GetNRows(sec); rowi++) {
+ for (UInt_t padi=0; padi<roc->GetNPads(sec,rowi); padi++) {
+ roc->GetPositionGlobal(sec,rowi,padi,xyz);
+ binx = 1+TMath::Nint((xyz[0]+250.)*0.5);
+ biny = 1+TMath::Nint((xyz[1]+250.)*0.5);
+ his2D->SetBinContent(binx,biny,npulse);
+ }
+ }
+ } else {
+ roc->GetPositionGlobal(sec,row,pad,xyz);
+ binx = 1+TMath::Nint((xyz[0]+250.)*0.5);
+ biny = 1+TMath::Nint((xyz[1]+250.)*0.5);
+
+ his2D->SetBinContent(binx,biny,npulse);
+ }
+ if (iHist%100==0){ printf("hist %d out of %d\n",iHist,nHist);}
+ }
+ his2D->SetXTitle("x (cm)");
+ his2D->SetYTitle("y (cm)");
+ his2D->SetStats(0);
+
+ his2D->DrawCopy("colz");
+
+ if (!side) {
+ gPad->SetTitle("A side");
+ } else {
+ gPad->SetTitle("C side");
+ }
+
+ return his2D;
+}
+
//____________________________________________________________________________
-TNtuple *AliTPCCalibTCF::PlotOccupSummary(const char *nameFile, Int_t nPulseMin) {
+void AliTPCCalibTCF::PlotOccupSummary(const char *nameFile, Int_t side, Int_t nPulseMin) {
//
// Plots the number of summed pulses per pad above a given minimum at the
- // pad position
+ // pad position at a given TPC side
// 'nameFile': root-file created with the Process function
//
TFile *file = new TFile(nameFile,"READ");
-
TH1F *his;
TKey *key;
TIter next( file->GetListOfKeys() );
+
+ char nameFileOut[100];
+ snprintf(nameFileOut,100,"Occup-%s",nameFile);
+ TFile fileOut(nameFileOut,"RECREATE");
+ // fileOut.cd();
+
TNtuple *ntuple = new TNtuple("ntuple","ntuple","x:y:z:npulse");
+ // ntuple->SetDirectory(0); // force to be memory resistent
+
+ Int_t nHist=file->GetNkeys();
+ if (!nHist) { return; }
+ Int_t iHist = 0;
+
+ Int_t secWise = 0;
- Int_t nPads = 0;
while ((key = (TKey *) next())) { // loop over histograms within the file
+
+ TString name(key->GetName());
+ if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl
his = (TH1F*)file->Get(key->GetName()); // copy object to memory
-
+ iHist++;
Int_t npulse = (Int_t)his->GetBinContent(1);
Int_t sec = (Int_t)his->GetBinContent(2);
Int_t row = (Int_t)his->GetBinContent(3);
Int_t pad = (Int_t)his->GetBinContent(4);
- if (row==-1 & pad==-1) { // summed pulses per sector
+ if ( (row<0) && (pad<0) ) { // row and pad are equal to -1, then -> summed pulses per sector
row = 40; pad = 40; // set to approx middle row for better plot
+ secWise=1;
}
Float_t *pos = new Float_t[3];
AliTPCROC::Instance()->GetPositionGlobal(sec,row,pad,pos);
if (npulse>=nPulseMin) {
ntuple->Fill(pos[0],pos[1],pos[2],npulse);
- printf("%d collected pulses in sector %d row %d pad %d\n",npulse,sec,row,pad);
+ if (iHist%100==0){ printf("hist %d out of %d\n",iHist,nHist);}
}
pos->~Float_t();
- nPads++;
}
-
- TCanvas *c1 = new TCanvas("TCanvas","Number of pulses found",1000,500);
- c1->Divide(2,1);
- char cSel[100];
- gStyle->SetPalette(1);
- gStyle->SetLabelOffset(-0.03,"Z");
- if (nPads<72) { // pulse per pad
+ if (secWise) { // pulse per sector
ntuple->SetMarkerStyle(8);
ntuple->SetMarkerSize(4);
- } else { // pulse per sector
+ } else { // pulse per Pad
ntuple->SetMarkerStyle(7);
}
- c1->cd(1);
- sprintf(cSel,"z>0&&npulse>=%d",nPulseMin);
- ntuple->Draw("y:x:npulse",cSel,"colz");
- gPad->SetTitle("A side");
+ char cSel[100];
+ if (!side) {
+ snprintf(cSel,100,"z>0&&npulse>=%d",nPulseMin);
+ ntuple->Draw("y:x:npulse",cSel,"colz");
+ } else {
+ snprintf(cSel,100,"z<0&&npulse>=%d",nPulseMin);
+ ntuple->Draw("y:x:npulse",cSel,"colz");
+ }
- c1->cd(2);
- sprintf(cSel,"z<0&&npulse>%d",nPulseMin);
- ntuple->Draw("y:x:npulse",cSel,"colz");
- gPad->SetTitle("C side");
+ if (!side) {
+ gPad->SetTitle("A side");
+ } else {
+ gPad->SetTitle("C side");
+ }
- file->Close();
- return ntuple;
+ ntuple->Write();
+ fileOut.Close();
+ file->Close();
}
//____________________________________________________________________________
TFile file(nameFileQuality,"READ");
TNtuple *qualityTuple = (TNtuple*)file.Get("TCFquality");
- gStyle->SetPalette(1);
- qualityTuple->Draw(plotSpec,cut,pOpt);
+ //gStyle->SetPalette(1);
-}
-
-//____________________________________________________________________________
-void AliTPCCalibTCF::DumpTCFparamToFile(const char *nameFileTCF,const char *nameFileOut)
-{
- //
- // Writes the TCF parameters from file 'nameFileTCF' to a output file
- //
-
- // Note: currently just TCF parameters per Sector or TCF parameters for pad
- // which were analyzed. There is no method included so far to export
- // parameters for not analyzed pad, which means there are eventually
- // missing TCF parameters
- // TODO: carefull! Fill up missing pads with averaged (sector) values?
+ TH2F *his2D = new TH2F(plotSpec,nameFileQuality,11,-10,1,25,1,100);
+ char plSpec[100];
+ snprintf(plSpec,100,"%s>>%s",plotSpec,plotSpec);
+ qualityTuple->Draw(plSpec,cut,pOpt);
+
+ gStyle->SetLabelSize(0.03,"X");
+ gStyle->SetLabelSize(0.03,"Y");
+ gStyle->SetLabelSize(0.03,"Z");
+ gStyle->SetLabelOffset(-0.02,"X");
+ gStyle->SetLabelOffset(-0.01,"Y");
+ gStyle->SetLabelOffset(-0.03,"Z");
+ his2D->GetXaxis()->SetTitle("max. undershot [ADC]");
+ his2D->GetYaxis()->SetTitle("width Reduction [%]");
- // open file with TCF parameters
- TFile fileTCF(nameFileTCF,"READ");
- TNtuple *paramTuple = (TNtuple*)fileTCF.Get("TCFparam");
-
- // open output txt file ...
- FILE *output;
- output=fopen(nameFileOut,"w"); // open outfile.
+ his2D->DrawCopy(pOpt);
- // Header line
- Int_t sectorWise = paramTuple->GetEntries("row==-1&&pad==-1");
- if (sectorWise) {
- fprintf(output,"sector \t Z0 \t\t Z1 \t\t Z2 \t\t P0 \t\t P1 \t\t P2\n");
- } else {
- fprintf(output,"sector \t row \t pad \t Z0 \t\t Z1 \t\t Z2 \t\t P0 \t\t P1 \t\t P2\n");
- }
+ gPad->SetPhi(0.1);gPad->SetTheta(90);
- for (Int_t i=0; i<paramTuple->GetEntries(); i++) {
- paramTuple->GetEntry(i);
- Float_t *p = paramTuple->GetArgs();
-
- // _______________________________________________________________
- // to Tuple to txt file - unsorted printout
-
- for (Int_t i=0; i<10; i++){
- if (sectorWise) {
- if (i<1) fprintf(output,"%3.0f \t ",p[i]); // sector info
- if (i>3) fprintf(output,"%1.4f \t ",p[i]); // TCF param
- } else {
- if (i<3) fprintf(output,"%3.0f \t ",p[i]); // pad info
- if (i>3) fprintf(output,"%1.4f \t ",p[i]); // TCF param
- }
- }
- fprintf(output,"\n");
- }
-
- // close output txt file
- fprintf(output,"\n");
- fclose(output);
+ his2D->~TH2F();
- fileTCF.Close();
-
-
}
-
-
//_____________________________________________________________________________
Int_t AliTPCCalibTCF::FitPulse(TNtuple *dataTuple, Double_t *coefZ, Double_t *coefP) {
//
// initialize TMinuit with a maximum of 8 params
- TMinuit *gMinuit = new TMinuit(8);
- gMinuit->mncler(); // Reset Minuit's list of paramters
- gMinuit->SetPrintLevel(-1); // No Printout
- gMinuit->SetFCN(AliTPCCalibTCF::FitFcn); // To set the address of the
+ TMinuit *minuitFit = new TMinuit(8);
+ minuitFit->mncler(); // Reset Minuit's list of paramters
+ minuitFit->SetPrintLevel(-1); // No Printout
+ minuitFit->SetFCN(AliTPCCalibTCF::FitFcn); // To set the address of the
// minimization function
- gMinuit->SetObjectFit(dataTuple);
+ minuitFit->SetObjectFit(dataTuple);
Double_t arglist[10];
Int_t ierflg = 0;
arglist[0] = 1;
- gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
+ minuitFit->mnexcm("SET ERR", arglist ,1,ierflg);
// Set standard starting values and step sizes for each parameter
// upper and lower limit (in a reasonable range) are set to improve
static Double_t min[8] = {100, 3., 0.1, 0.2, 3., 60., 0., 2.0};
static Double_t max[8] = {200, 20., 5., 3., 30., 300., 20., 2.5};
- gMinuit->mnparm(0, "A1", vstart[0], step[0], min[0], max[0], ierflg);
- gMinuit->mnparm(1, "A2", vstart[1], step[1], min[1], max[1], ierflg);
- gMinuit->mnparm(2, "A3", vstart[2], step[2], min[2], max[2], ierflg);
- gMinuit->mnparm(3, "T1", vstart[3], step[3], min[3], max[3], ierflg);
- gMinuit->mnparm(4, "T2", vstart[4], step[4], min[4], max[4], ierflg);
- gMinuit->mnparm(5, "T3", vstart[5], step[5], min[5], max[5], ierflg);
- gMinuit->mnparm(6, "T0", vstart[6], step[6], min[6], max[6], ierflg);
- gMinuit->mnparm(7, "TTP", vstart[7], step[7], min[7], max[7],ierflg);
- gMinuit->FixParameter(7); // 2.24 ... out of pulserRun Fit (->IRF)
+ minuitFit->mnparm(0, "A1", vstart[0], step[0], min[0], max[0], ierflg);
+ minuitFit->mnparm(1, "A2", vstart[1], step[1], min[1], max[1], ierflg);
+ minuitFit->mnparm(2, "A3", vstart[2], step[2], min[2], max[2], ierflg);
+ minuitFit->mnparm(3, "T1", vstart[3], step[3], min[3], max[3], ierflg);
+ minuitFit->mnparm(4, "T2", vstart[4], step[4], min[4], max[4], ierflg);
+ minuitFit->mnparm(5, "T3", vstart[5], step[5], min[5], max[5], ierflg);
+ minuitFit->mnparm(6, "T0", vstart[6], step[6], min[6], max[6], ierflg);
+ minuitFit->mnparm(7, "TTP", vstart[7], step[7], min[7], max[7],ierflg);
+ minuitFit->FixParameter(7); // 2.24 ... out of pulserRun Fit (->IRF)
// Now ready for minimization step
arglist[0] = 2000; // max num of iterations
arglist[1] = 0.1; // tolerance
- gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
+ minuitFit->mnexcm("MIGRAD", arglist ,2,ierflg);
Double_t p1 = 0.0 ;
- gMinuit->mnexcm("SET NOW", &p1 , 0, ierflg) ; // No Warnings
+ minuitFit->mnexcm("SET NOW", &p1 , 0, ierflg) ; // No Warnings
if (ierflg == 4) { // Fit failed
for (Int_t i=0;i<3;i++) {
coefP[i] = 0;
coefZ[i] = 0;
}
- gMinuit->~TMinuit();
+ minuitFit->~TMinuit();
return 0;
} else { // Fit successfull
for (Int_t i=0;i<6;i++) {
Double_t err = 0;
Double_t val = 0;
- gMinuit->GetParameter(i,val,err);
+ minuitFit->GetParameter(i,val,err);
fitParam[i] = val;
}
// TCF coefficients which are used for the equalisation step (stage)
// ZERO/POLE Filter
- coefZ[0] = exp(-1/valuePZ[2]);
- coefZ[1] = exp(-1/valuePZ[3]);
- coefP[0] = exp(-1/valuePZ[0]);
- coefP[1] = exp(-1/valuePZ[1]);
+ coefZ[0] = TMath::Exp(-1/valuePZ[2]);
+ coefZ[1] = TMath::Exp(-1/valuePZ[3]);
+ coefP[0] = TMath::Exp(-1/valuePZ[0]);
+ coefP[1] = TMath::Exp(-1/valuePZ[1]);
fitParam->~Double_t();
valuePZ->~Double_t();
- gMinuit->~TMinuit();
+ minuitFit->~TMinuit();
return 1;
//____________________________________________________________________________
-void AliTPCCalibTCF::FitFcn(Int_t &/*nPar*/, Double_t */*grad*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
+void AliTPCCalibTCF::FitFcn(Int_t &/*nPar*/, Double_t */*grad*/, Double_t &f, Double_t * const par, Int_t /*iflag*/)
{
//
// Minimization function needed for TMinuit with FitFunction included
if (t<0) {
sigFit = 0;
} else {
- Double_t f1 = 1/pow((4-ttp/par[3]),5)*(24*ttp*exp(4)*(exp(-t/par[3]) - exp(-4*t/ttp) * ( 1+t*(4-ttp/par[3])/ttp+pow(t*(4-ttp/par[3])/ttp,2)/2 + pow(t*(4-ttp/par[3])/ttp,3)/6 + pow(t*(4-ttp/par[3])/ttp,4)/24)));
+ Double_t f1 = 1/TMath::Power((4-ttp/par[3]),5)*(24*ttp*TMath::Exp(4)*(TMath::Exp(-t/par[3]) - TMath::Exp(-4*t/ttp) * ( 1+t*(4-ttp/par[3])/ttp+TMath::Power(t*(4-ttp/par[3])/ttp,2)/2 + TMath::Power(t*(4-ttp/par[3])/ttp,3)/6 + TMath::Power(t*(4-ttp/par[3])/ttp,4)/24)));
- Double_t f2 = 1/pow((4-ttp/par[4]),5)*(24*ttp*exp(4)*(exp(-t/par[4]) - exp(-4*t/ttp) * ( 1+t*(4-ttp/par[4])/ttp+pow(t*(4-ttp/par[4])/ttp,2)/2 + pow(t*(4-ttp/par[4])/ttp,3)/6 + pow(t*(4-ttp/par[4])/ttp,4)/24)));
+ Double_t f2 = 1/TMath::Power((4-ttp/par[4]),5)*(24*ttp*TMath::Exp(4)*(TMath::Exp(-t/par[4]) - TMath::Exp(-4*t/ttp) * ( 1+t*(4-ttp/par[4])/ttp+TMath::Power(t*(4-ttp/par[4])/ttp,2)/2 + TMath::Power(t*(4-ttp/par[4])/ttp,3)/6 + TMath::Power(t*(4-ttp/par[4])/ttp,4)/24)));
- Double_t f3 = 1/pow((4-ttp/par[5]),5)*(24*ttp*exp(4)*(exp(-t/par[5]) - exp(-4*t/ttp) * ( 1+t*(4-ttp/par[5])/ttp+pow(t*(4-ttp/par[5])/ttp,2)/2 + pow(t*(4-ttp/par[5])/ttp,3)/6 + pow(t*(4-ttp/par[5])/ttp,4)/24)));
+ Double_t f3 = 1/TMath::Power((4-ttp/par[5]),5)*(24*ttp*TMath::Exp(4)*(TMath::Exp(-t/par[5]) - TMath::Exp(-4*t/ttp) * ( 1+t*(4-ttp/par[5])/ttp+TMath::Power(t*(4-ttp/par[5])/ttp,2)/2 + TMath::Power(t*(4-ttp/par[5])/ttp,3)/6 + TMath::Power(t*(4-ttp/par[5])/ttp,4)/24)));
sigFit = par[0]*f1 + par[1]*f2 +par[2]*f3;
}
// to the different stages of the TCF filter
// (e.g. first 2 fit parameters represent the electron signal itself!)
- if (param[3]==param[4]) {param[3]=param[3]+0.0001;}
- if (param[5]==param[4]) {param[5]=param[5]+0.0001;}
+ if ((param[3]-param[4]) <1e-5 ) {param[3]=param[3]+0.0001;} // if equal
+ if ((param[5]-param[4]) <1e-5 ) {param[5]=param[5]+0.0001;} // if equal
- if ((param[5]>param[4])&(param[5]>param[3])) {
+ if ((param[5]>param[4])&&(param[5]>param[3])) {
if (param[4]>=param[3]) {
vA1 = param[0]; vA2 = param[1]; vA3 = param[2];
vTT1 = param[3]; vTT2 = param[4]; vTT3 = param[5];
vA1 = param[1]; vA2 = param[0]; vA3 = param[2];
vTT1 = param[4]; vTT2 = param[3]; vTT3 = param[5];
}
- } else if ((param[4]>param[5])&(param[4]>param[3])) {
+ } else if ((param[4]>param[5])&&(param[4]>param[3])) {
if (param[5]>=param[3]) {
vA1 = param[0]; vA2 = param[2]; vA3 = param[1];
vTT1 = param[3]; vTT2 = param[5]; vTT3 = param[4];
vA1 = param[2]; vA2 = param[0]; vA3 = param[1];
vTT1 = param[5]; vTT2 = param[3]; vTT3 = param[4];
}
- } else if ((param[3]>param[4])&(param[3]>param[5])) {
+ } else if ((param[3]>param[4])&&(param[3]>param[5])) {
if (param[5]>=param[4]) {
vA1 = param[1]; vA2 = param[2]; vA3 = param[0];
vTT1 = param[4]; vTT2 = param[5]; vTT3 = param[3];
//____________________________________________________________________________
-void AliTPCCalibTCF::Equalization(TNtuple *dataTuple, Double_t *coefZ, Double_t *coefP) {
+Int_t AliTPCCalibTCF::Equalization(TNtuple *dataTuple, Double_t *coefZ, Double_t *coefP) {
//
// calculates the 3rd set of TCF parameters (remaining 2 PZ values) in
// order to restore the original pulse height and adds them to the passed arrays
//
- Double_t *s0 = new Double_t[1000]; // original pulse
- Double_t *s1 = new Double_t[1000]; // pulse after 1st PZ filter
- Double_t *s2 = new Double_t[1000]; // pulse after 2nd PZ filter
-
const Int_t kPulseLength = dataTuple->GetEntries();
-
+
+ if (kPulseLength<2) {
+ // prinft("PulseLength does not make sense\n");
+ return 0;
+ }
+
+ Double_t *s0 = new Double_t[kPulseLength]; // original pulse
+ Double_t *s1 = new Double_t[kPulseLength]; // pulse after 1st PZ filter
+ Double_t *s2 = new Double_t[kPulseLength]; // pulse after 2nd PZ filter
+
for (Int_t ipos=0; ipos<kPulseLength; ipos++) {
dataTuple->GetEntry(ipos);
Float_t *p = dataTuple->GetArgs();
coefP[2] = 0;
coefZ[2] = (s2ampl - s0ampl)/s0[s0pos-1];
} else { // same height ? will most likely not happen ?
+ printf("No equalization because of identical height\n");
coefP[2] = 0;
coefZ[2] = 0;
}
- s0->~Double_t();
- s1->~Double_t();
- s2->~Double_t();
-
+ delete [] s0;
+ delete [] s1;
+ delete [] s2;
+
+ // if equalization out of range (<0 or >=1) it failed!
+ // if ratio of amplitudes of fittet to original pulse < 0.9 it failed!
+ if (coefP[2]<0 || coefZ[2]<0 || coefP[2]>=1 || coefZ[2]>=1 || TMath::Abs(s2ampl / s0ampl)<0.9) {
+ return 0;
+ } else {
+ return 1;
+ }
+
}
//____________________________________________________________________________
-Int_t AliTPCCalibTCF::FindCorTCFparam(TH1F *hisIn, const char *nameFileTCF, Double_t *coefZ, Double_t *coefP) {
+Int_t AliTPCCalibTCF::FindCorTCFparam(TH1F * const hisIn, const char *nameFileTCF, Double_t *coefZ, Double_t *coefP) {
//
// This function searches for the correct TCF parameters to the given
// histogram 'hisIn' within the file 'nameFileTCF'
char sel[100];
if ( paramTuple->GetEntries("row==-1&&pad==-1") ) {
// parameters per SECTOR
- sprintf(sel,"sec==%d&&row==-1&&pad==-1",sector);
+ snprintf(sel,100,"sec==%d&&row==-1&&pad==-1",sector);
} else {
// parameters per PAD
- sprintf(sel,"sec==%d&&row==%d&&pad==%d",sector,row,pad);
+ snprintf(sel,100,"sec==%d&&row==%d&&pad==%d",sector,row,pad);
}
// list should contain just ONE entry! ... otherwise there is a mistake!
paramTuple->GetEntry(pos); // get specific TCF parameters
Float_t *p = paramTuple->GetArgs();
// check ...
- if(sector==p[0]) {printf("sector ok ... "); }
- if(row==p[1]) {printf("row ok ... "); }
- if(pad==p[2]) {printf("pad ok ... \n"); }
+ if((sector-p[0])<1e-5) {printf("sector ok ... "); }
+ if((row-p[1])<1e-5) {printf("row ok ... "); }
+ if((pad-p[2])<1e-5) {printf("pad ok ... \n"); }
// number of averaged pulses used to produce TCF params
nPulse = (Int_t)p[3];
} else { // no specific TCF parameters found for this pad
- printf("no specific TCF paramaters found for pad in ...\n");
- printf("in Sector %d | Row %d | Pad %d |\n", sector, row, pad);
+ printf(" no specific TCF paramaters found for pad in ...\n");
+ printf(" Sector %d | Row %d | Pad %d |\n", sector, row, pad);
nPulse = 0;
coefZ[0] = 0; coefP[0] = 0;
coefZ[1] = 0; coefP[1] = 0;
startOfPulse = 1;
posOfStart = i;
}
- if( (sig < threshADC) && (startOfPulse == 1) ){
+ if( (sig <= threshADC) && (startOfPulse == 1) ){
widthFound = 1;
width = i - posOfStart + 1;
}
startOfPulseTCF = 1;
posOfStartTCF = i;
}
- if( (sigTCF < threshADC) && (startOfPulseTCF == 1) ){
+ if( (sigTCF <= threshADC) && (startOfPulseTCF == 1) ){
widthFoundTCF = 1;
widthTCF = i -posOfStartTCF + 1;
}
}
// Search for maximal undershot (is equal to minimum after the pulse)
- if ( (undershotStart==1)&&(i<(posOfStartTCF+widthTCF+20)) ) {
+ if ( ((undershotStart-1)<1e-7)&&(i<(posOfStartTCF+widthTCF+20)) ) {
if (maxUndershot>sigTCF) { maxUndershot = sigTCF; }
}
for (Int_t ipos = 0; ipos<6; ipos++) {
// before the pulse
tempRMSHis->Fill(hisIn->GetBinContent(ipos+5));
- // at the end
+ // at the end
tempRMSHis->Fill(hisIn->GetBinContent(hisIn->GetNbinsX()-ipos));
}
Double_t pulseRMS = tempRMSHis->GetRMS();
printf("width reduction [percent]:\t\t %3.1f\n", widthReduct);
printf("mean undershot [ADC]:\t\t\t %3.1f\n", undershot);
printf("maximum of undershot after pulse [ADC]: %3.1f\n", maxUndershot);
- printf("RMS of the original pulse [ADC]: \t %3.2f\n\n", pulseRMS);
+ printf("RMS of the original (or summed) pulse [ADC]: \t %3.2f\n\n", pulseRMS);
}
//____________________________________________________________________________
-TNtuple *AliTPCCalibTCF::ApplyTCFilter(TH1F *hisIn, Double_t *coefZ, Double_t *coefP, Int_t plotFlag) {
+TNtuple *AliTPCCalibTCF::ApplyTCFilter(TH1F * const hisIn, Double_t * const coefZ, Double_t * const coefP, Int_t plotFlag) {
//
// Applies the given TCF parameters on the given pulse via the ALTRO emulator
// class (discret values) and stores both pulses into a returned TNtuple
}
// transform TCF parameters into ALTRO readable format (Integer)
- Int_t* valP = new Int_t[3];
- Int_t* valZ = new Int_t[3];
+ Int_t valK[3];
+ Int_t valL[3];
for (Int_t i=0; i<3; i++) {
- valP[i] = (Int_t)(coefP[i]*(pow(2,16)-1));
- valZ[i] = (Int_t)(coefZ[i]*(pow(2,16)-1));
+ valK[i] = (Int_t)(coefP[i]*(TMath::Power(2,16)-1));
+ valL[i] = (Int_t)(coefZ[i]*(TMath::Power(2,16)-1));
}
// discret ALTRO EMULATOR ____________________________
AliTPCAltroEmulator *altro = new AliTPCAltroEmulator(nbins, signalOutD);
altro->ConfigAltro(0,1,0,0,0,0); // perform just the TailCancelation
- altro->ConfigTailCancellationFilter(valP[0],valP[1],valP[2],valZ[0],valZ[1],valZ[2]);
+ altro->ConfigTailCancellationFilter(valK[0],valK[1],valK[2],valL[0],valL[1],valL[2]);
altro->RunEmulation();
delete altro;
if (plotFlag) {
char hname[100];
- sprintf(hname,"sec%drow%dpad%d",sector,row,pad);
+ snprintf(hname,100,"sec%drow%dpad%d",sector,row,pad);
new TCanvas(hname,hname,600,400);
//just plotting non-discret pulses | they look pretties in case of mean sig ;-)
pulseTuple->Draw("sigND:timebin","","L");
// pulseTuple->Draw("sigAfterTCF:timebin","","Lsame");
}
- valP->~Int_t();
- valZ->~Int_t();
-
- signalIn->~Double_t();
- signalOut->~Double_t();
- delete signalIn;
- delete signalOut;
-
+ delete [] signalIn;
+ delete [] signalOut;
+ delete [] signalInD;
+ delete [] signalOutD;
+
return pulseTuple;
}
-
-
//____________________________________________________________________________
void AliTPCCalibTCF::PrintPulseThresholds() {
//
printf(" %4.0d [ADC] ... lower pulse height limit \n", fLowPulseLim);
printf(" %4.0d [ADC] ... upper pulse height limit \n", fUpPulseLim);
printf(" %4.1f [ADC] ... maximal pulse RMS \n", fRMSLim);
+ printf(" %4.1f [ADC] ... pulse/tail integral ratio \n", fRatioIntLim);
}
+
+
+//____________________________________________________________________________
+void AliTPCCalibTCF::MergeHistoPerFile(const char *fileNameIn, const char *fileNameSum, Int_t mode)
+{
+ // Gets histograms from fileNameIn and adds contents to fileSum
+ //
+ // If fileSum doesn't exist, fileSum is created
+ // mode = 0, just ONE BIG FILE ('fileSum') will be used
+ // mode = 1, one file per sector ('fileSum-Sec#.root') will be used
+ // mode=1 is much faster, but the additional function 'MergeToOneFile' has to be used in order to
+ // get one big and complete collection file again ...
+ //
+ // !Make sure not to add the same file more than once!
+
+ TFile fileIn(fileNameIn,"READ");
+ TH1F *hisIn;
+ TKey *key;
+ TIter next(fileIn.GetListOfKeys());
+ // opens a file, although, it might not be uses (see "mode")
+ TFile *fileOut = new TFile(fileNameSum,"UPDATE");
+ //fileOut.cd();
+
+ Int_t nHist=fileIn.GetNkeys();
+ Int_t iHist=0;
+
+ Int_t secPrev = -1;
+ char fileNameSumSec[100];
+
+
+ while((key=(TKey*)next())) {
+ const char *hisName = key->GetName();
+
+ TString name(key->GetName());
+ if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl
+
+ hisIn=(TH1F*)fileIn.Get(hisName);
+
+ Int_t numPulse=(Int_t)hisIn->GetBinContent(1);
+ Int_t sec=(Int_t)hisIn->GetBinContent(2);
+ Int_t pulseLength= hisIn->GetNbinsX()-4;
+
+ // in case of mode 1, store histos in files per sector
+ if (sec!=secPrev && mode != 0) {
+ if (secPrev>0) { // closing old file
+ fileOut->Close();
+ }
+ // opening new file
+ snprintf(fileNameSumSec,100,"%s-Sec%d.root",fileNameSum,sec);
+ fileOut = new TFile(fileNameSumSec,"UPDATE");
+ secPrev = sec;
+ }
+
+ // search for existing histogram
+ TH1F *his=(TH1F*)fileOut->Get(hisName);
+ if (iHist%100==0) {
+ printf("Histogram %d / %d, %s, Action: ",iHist,nHist,hisName);
+ if (!his) {
+ printf("NEW\n");
+ } else {
+ printf("ADD\n");
+ }
+ }
+ iHist++;
+
+ if (!his) {
+ his=hisIn;
+ his->Write(hisName);
+ } else {
+ his->AddBinContent(1,numPulse);
+ for (Int_t ii=5; ii<pulseLength+5; ii++) {
+ his->AddBinContent(ii,hisIn->GetBinContent(ii));
+ }
+ his->Write(hisName,TObject::kOverwrite);
+ }
+ }
+
+ printf("closing files (may take a while)...\n");
+ fileOut->Close();
+
+
+ fileIn.Close();
+ printf("...DONE\n\n");
+}
+
+
+//____________________________________________________________________________
+void AliTPCCalibTCF::MergeToOneFile(const char *nameFileSum) {
+
+ // Merges all Sec-files together ...
+ // this is an additional functionality for the function MergeHistsPerFile
+ // if for example mode=1
+
+ TH1F *hisIn;
+ TKey *key;
+
+ // just delete the file entries ...
+ TFile fileSumD(nameFileSum,"RECREATE");
+ fileSumD.Close();
+
+ char nameFileSumSec[100];
+
+ for (Int_t sec=0; sec<72; sec++) { // loop over all possible filenames
+
+ snprintf(nameFileSumSec,100,"%s-Sec%d.root",nameFileSum,sec);
+ TFile *fileSumSec = new TFile(nameFileSumSec,"READ");
+
+ Int_t nHist=fileSumSec->GetNkeys();
+ Int_t iHist=0;
+
+ if (nHist) { // file found \ NKeys not empty
+
+ TFile fileSum(nameFileSum,"UPDATE");
+ fileSum.cd();
+
+ printf("Sector file %s found\n",nameFileSumSec);
+ TIter next(fileSumSec->GetListOfKeys());
+ while( (key=(TKey*)next()) ) {
+ const char *hisName = key->GetName();
+ TString name(hisName);
+ if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl
+ hisIn=(TH1F*)fileSumSec->Get(hisName);
+
+
+ if (iHist%100==0) {
+ printf("found histogram %d / %d, %s\n",iHist,nHist,hisName);
+ }
+ iHist++;
+
+ // TH1F *his = (TH1F*)hisIn->Clone(hisName);
+ hisIn->Write(hisName);
+
+ }
+ printf("Saving histograms from sector %d (may take a while) ...",sec);
+ fileSum.Close();
+
+ }
+ fileSumSec->Close();
+ }
+ printf("...DONE\n\n");
+}
+
+
+//____________________________________________________________________________
+Int_t AliTPCCalibTCF::DumpTCFparamToFilePerPad(const char *nameFileTCFPerPad,const char *nameFileTCFPerSec, const char *nameMappingFile) {
+ //
+ // Writes TCF parameters per PAD to .data file
+ //
+ // from now on: "roc" refers to the offline sector numbering
+ // "sector" refers to the 18 sectors per side
+ //
+ // Gets TCF parameters of single pads from nameFileTCFPerPad and writes them to
+ // the file 'tpcTCFparamPAD.data'
+ //
+ // If there are parameters for a pad missing, then the parameters of the roc,
+ // in which the pad is located, are used as the pad parameters. The parameters for
+ // the roc are retreived from nameFileTCFPerSec. If there are parameters for
+ // a roc missing, then the parameters are set to -1.
+
+ Float_t k0 = -1, k1 = -1, k2 = -1, l0 = -1, l1 = -1, l2 = -1;
+ Int_t roc, row, pad, side, sector, rcu, hwAddr;
+ Int_t entryNum = 0;
+ Int_t checksum = 0;
+ Int_t tpcPadNum = 557568;
+ Int_t validFlag = 1; // 1 if parameters for pad exist, 0 if they are only inherited from the roc
+
+ // get file/tuple with parameters per pad
+ TFile fileTCFparam(nameFileTCFPerPad);
+ TNtuple *paramTuple = (TNtuple*)fileTCFparam.Get("TCFparam");
+
+ // get mapping file
+ // usual location of mapping file: $ALICE_ROOT/TPC/Calib/tpcMapping.root
+ TFile *fileMapping = new TFile(nameMappingFile, "read");
+ AliTPCmapper *mapping = (AliTPCmapper*) fileMapping->Get("tpcMapping");
+ delete fileMapping;
+
+ if (mapping == 0) {
+ printf("Failed to get mapping object from %s. ...\n", nameMappingFile);
+ return -1;
+ } else {
+ printf("Got mapping object from %s\n", nameMappingFile);
+ }
+
+ Bool_t *entryID = new Bool_t[7200000]; // helping vector
+ for (Int_t ii = 0; ii<7200000; ii++) {
+ entryID[ii]=0;
+ }
+
+ // creating outputfile
+ ofstream fileOut;
+ char nameFileOut[255];
+ snprintf(nameFileOut,255,"tpcTCFparamPAD.data");
+ fileOut.open(nameFileOut);
+ // following not used:
+ // char headerLine[255];
+ // snprintf(headerLine,255,"15\tside\tsector\tRCU\tHWadr\tk0\tk1\tk2\tl0\tl1\tl2\tValidFlag");
+ // fileOut << headerLine << std::endl;
+ fileOut << "15" << std::endl;
+
+ // loop over nameFileTCFPerPad, write parameters into outputfile
+ // NOTE: NO SPECIFIC ORDER !!!
+ printf("\nstart assigning parameters to pad...\n");
+ for (Int_t iParam = 0; iParam < paramTuple->GetEntries(); iParam++) {
+ paramTuple->GetEntry(iParam);
+ Float_t *paramArgs = paramTuple->GetArgs();
+ roc = Int_t(paramArgs[0]);
+ row = Int_t(paramArgs[1]);
+ pad = Int_t(paramArgs[2]);
+ side = Int_t(mapping->GetSideFromRoc(roc));
+ sector = Int_t(mapping->GetSectorFromRoc(roc));
+ rcu = Int_t(mapping->GetRcu(roc,row,pad));
+ hwAddr = Int_t(mapping->GetHWAddress(roc,row,pad));
+ k0 = TMath::Nint(paramArgs[7] * (TMath::Power(2,16) - 1));
+ k1 = TMath::Nint(paramArgs[8] * (TMath::Power(2,16) - 1));
+ k2 = TMath::Nint(paramArgs[9] * (TMath::Power(2,16) - 1));
+ l0 = TMath::Nint(paramArgs[4] * (TMath::Power(2,16) - 1));
+ l1 = TMath::Nint(paramArgs[5] * (TMath::Power(2,16) - 1));
+ l2 = TMath::Nint(paramArgs[6] * (TMath::Power(2,16) - 1));
+ if (entryNum%10000==0) {
+ printf("assigned pad %i / %i\n",entryNum,tpcPadNum);
+ }
+
+ fileOut << entryNum++ << "\t" << side << "\t" << sector << "\t" << rcu << "\t" << hwAddr << "\t";
+ fileOut << k0 << "\t" << k1 << "\t" << k2 << "\t" << l0 << "\t" << l1 << "\t" << l2 << "\t" << validFlag << std::endl;
+ entryID[roc*100000 + row*1000 + pad] = 1;
+ }
+
+ // Wrote all found TCF params per pad into data file
+ // NOW FILLING UP THE REST WITH THE PARAMETERS FROM THE ROC MEAN
+
+ // get file/tuple with parameters per roc
+ TFile fileSecTCFparam(nameFileTCFPerSec);
+ TNtuple *paramTupleSec = (TNtuple*)fileSecTCFparam.Get("TCFparam");
+
+ // loop over all pads and get/write parameters for pads which don't have
+ // parameters assigned yet
+ validFlag = 0;
+ for (roc = 0; roc<72; roc++) {
+ side = Int_t(mapping->GetSideFromRoc(roc));
+ sector = Int_t(mapping->GetSectorFromRoc(roc));
+ for (Int_t iParamSec = 0; iParamSec < paramTupleSec->GetEntries(); iParamSec++) {
+ paramTupleSec->GetEntry(iParamSec);
+ Float_t *paramArgsSec = paramTupleSec->GetArgs();
+ if ((paramArgsSec[0]-roc)<1e-7) { // if roc is found
+ k0 = TMath::Nint(paramArgsSec[7] * (TMath::Power(2,16) - 1));
+ k1 = TMath::Nint(paramArgsSec[8] * (TMath::Power(2,16) - 1));
+ k2 = TMath::Nint(paramArgsSec[9] * (TMath::Power(2,16) - 1));
+ l0 = TMath::Nint(paramArgsSec[4] * (TMath::Power(2,16) - 1));
+ l1 = TMath::Nint(paramArgsSec[5] * (TMath::Power(2,16) - 1));
+ l2 = TMath::Nint(paramArgsSec[6] * (TMath::Power(2,16) - 1));
+ break;
+ } else {
+ k0 = k1 = k2 = l0 = l1 = l2 = -1;
+ }
+ }
+ for (row = 0; row<mapping->GetNpadrows(roc); row++) {
+ for (pad = 0; pad<mapping->GetNpads(roc,row); pad++) {
+ if (entryID[roc*100000 + row*1000 + pad]==1) {
+ continue;
+ }
+
+ entryID[roc*100000 + row*1000 + pad] = 1;
+ rcu = Int_t(mapping->GetRcu(roc,row,pad));
+ hwAddr = Int_t(mapping->GetHWAddress(roc,row,pad));
+ if (entryNum%10000==0) {
+ printf("assigned pad %i / %i\n",entryNum,tpcPadNum);
+ }
+
+ fileOut << entryNum++ << "\t" << side << "\t" << sector << "\t" << rcu << "\t" << hwAddr << "\t";
+ fileOut << k0 << "\t" << k1 << "\t" << k2 << "\t" << l0 << "\t" << l1 << "\t" << l2 << "\t" << validFlag << std::endl;
+ }
+ }
+ }
+
+ printf("assigned pad %i / %i\ndone assigning\n",entryNum,tpcPadNum);
+
+ // check if correct amount of sets of parameters were written
+ for (Int_t ii = 0; ii<7200000; ii++) {
+ checksum += entryID[ii];
+ }
+ if (checksum == tpcPadNum) {
+ printf("checksum ok, sets of parameters written = %i\n",checksum);
+ } else {
+ printf("\nCHECKSUM WRONG, sets of parameters written = %i, should be %i\n\n",checksum,tpcPadNum);
+ }
+
+ // closing & destroying
+ fileOut.close();
+ fileTCFparam.Close();
+ fileSecTCFparam.Close();
+ delete [] entryID;
+ printf("output written to file: %s\n",nameFileOut);
+ return 0;
+}
+
+
+
+//____________________________________________________________________________
+Int_t AliTPCCalibTCF::DumpTCFparamToFilePerSector(const char *nameFileTCFPerSec, const char *nameMappingFile) {
+ //
+ // Writes TCF parameters per SECTOR (=ROC) to .data file
+ //
+ // from now on: "roc" refers to the offline sector numbering
+ // "sector" refers to the 18 sectors per side
+ //
+ // Gets TCF parameters of a roc from nameFileTCFPerSec and writes them to
+ // the file 'tpcTCFparamSector.data'
+ //
+ // If there are parameters for a roc missing, then the parameters are set to -1
+
+ Float_t k0 = -1, k1 = -1, k2 = -1, l0 = -1, l1 = -1, l2 = -1;
+ Int_t entryNum = 0;
+ Int_t validFlag = 0; // 1 if parameters for roc exist
+
+ // get file/tuple with parameters per roc
+ TFile fileTCFparam(nameFileTCFPerSec);
+ TNtuple *paramTupleSec = (TNtuple*)fileTCFparam.Get("TCFparam");
+
+
+ // get mapping file
+ // usual location of mapping file: $ALICE_ROOT/TPC/Calib/tpcMapping.root
+ TFile *fileMapping = new TFile(nameMappingFile, "read");
+ AliTPCmapper *mapping = (AliTPCmapper*) fileMapping->Get("tpcMapping");
+ delete fileMapping;
+
+ if (mapping == 0) {
+ printf("Failed to get mapping object from %s. ...\n", nameMappingFile);
+ return -1;
+ } else {
+ printf("Got mapping object from %s\n", nameMappingFile);
+ }
+
+
+ // creating outputfile
+
+ ofstream fileOut;
+ char nameFileOut[255];
+ snprintf(nameFileOut,255,"tpcTCFparamSector.data");
+ fileOut.open(nameFileOut);
+ // following not used:
+ // char headerLine[255];
+ // snprintf(headerLine,255,"16\tside\tsector\tRCU\tHWadr\tk0\tk1\tk2\tl0\tl1\tl2\tValidFlag");
+ // fileOut << headerLine << std::endl;
+ fileOut << "16" << std::endl;
+
+ // loop over all rcu's in the TPC (6 per sector)
+ printf("\nstart assigning parameters to rcu's...\n");
+ for (Int_t side = 0; side<2; side++) {
+ for (Int_t sector = 0; sector<18; sector++) {
+ for (Int_t rcu = 0; rcu<6; rcu++) {
+
+ validFlag = 0;
+ Int_t roc = Int_t(mapping->GetRocFromPatch(side, sector, rcu));
+
+ // get parameters (through loop search) for rcu from corresponding roc
+ for (Int_t iParam = 0; iParam < paramTupleSec->GetEntries(); iParam++) {
+ paramTupleSec->GetEntry(iParam);
+ Float_t *paramArgs = paramTupleSec->GetArgs();
+ if ((paramArgs[0]-roc)<1e-7) { // if roc is found
+ validFlag = 1;
+ k0 = TMath::Nint(paramArgs[7] * (TMath::Power(2,16) - 1));
+ k1 = TMath::Nint(paramArgs[8] * (TMath::Power(2,16) - 1));
+ k2 = TMath::Nint(paramArgs[9] * (TMath::Power(2,16) - 1));
+ l0 = TMath::Nint(paramArgs[4] * (TMath::Power(2,16) - 1));
+ l1 = TMath::Nint(paramArgs[5] * (TMath::Power(2,16) - 1));
+ l2 = TMath::Nint(paramArgs[6] * (TMath::Power(2,16) - 1));
+ break;
+ }
+ }
+ if (!validFlag) { // No TCF parameters found for this roc
+ k0 = k1 = k2 = l0 = l1 = l2 = -1;
+ }
+
+ fileOut << entryNum++ << "\t" << side << "\t" << sector << "\t" << rcu << "\t" << -1 << "\t";
+ fileOut << k0 << "\t" << k1 << "\t" << k2 << "\t" << l0 << "\t" << l1 << "\t" << l2 << "\t" << validFlag << std::endl;
+ }
+ }
+ }
+
+ printf("done assigning\n");
+
+ // closing files
+ fileOut.close();
+ fileTCFparam.Close();
+ printf("output written to file: %s\n",nameFileOut);
+ return 0;
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff