X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCCalibTCF.cxx;h=ba909d532324fd7b3597f8ff15e93836c56f3071;hb=dce050573d0e990f92bda56fc2f1a4fc7059ecc9;hp=be2f734df6a4bedfbbe34dfb31e9e2c6ad481e83;hpb=2c632057b20a755472de41d084d02d410a4d451d;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCCalibTCF.cxx b/TPC/AliTPCCalibTCF.cxx index be2f734df6a..ba909d53232 100644 --- a/TPC/AliTPCCalibTCF.cxx +++ b/TPC/AliTPCCalibTCF.cxx @@ -19,7 +19,7 @@ // 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 // // // /////////////////////////////////////////////////////////////////////////////// @@ -32,27 +32,35 @@ #include #include #include +#include +#include #include #include #include - #include "AliRawReaderRoot.h" +#include "AliRawHLTManager.h" #include "AliTPCRawStream.h" +#include "AliTPCRawStreamV3.h" #include "AliTPCROC.h" #include "AliTPCAltroEmulator.h" +#include "AliTPCmapper.h" +#include + 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 @@ -60,14 +68,15 @@ AliTPCCalibTCF::AliTPCCalibTCF() : } //_____________________________________________________________________________ -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 @@ -82,7 +91,8 @@ AliTPCCalibTCF::AliTPCCalibTCF(const AliTPCCalibTCF &tcf) : fPulseLength(tcf.fPulseLength), fLowPulseLim(tcf.fLowPulseLim), fUpPulseLim(tcf.fUpPulseLim), - fRMSLim(tcf.fRMSLim) + fRMSLim(tcf.fRMSLim), + fRatioIntLim(tcf.fRatioIntLim) { // // AliTPCCalibTCF copy constructor @@ -112,29 +122,269 @@ AliTPCCalibTCF::~AliTPCCalibTCF() // } + +//_____________________________________________________________________________ +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 (blGetSector(); + 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;iSetBinContent(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 && max10&& (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; iposGetBinContent(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; iposGetBinContent(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; iposGetBinContent(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; iposGetBinContent(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) { // @@ -146,10 +396,15 @@ void AliTPCCalibTCF::ProcessRawEvent(AliTPCRawStream *rawStream, const char *nam // 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(); @@ -157,7 +412,7 @@ void AliTPCCalibTCF::ProcessRawEvent(AliTPCRawStream *rawStream, const char *nam 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()) { @@ -165,31 +420,55 @@ void AliTPCCalibTCF::ProcessRawEvent(AliTPCRawStream *rawStream, const char *nam 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 @@ -197,9 +476,12 @@ void AliTPCCalibTCF::ProcessRawEvent(AliTPCRawStream *rawStream, const char *nam 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(); @@ -207,23 +489,43 @@ void AliTPCCalibTCF::ProcessRawEvent(AliTPCRawStream *rawStream, const char *nam 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 && maxlowLim && max= 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; iposGetBinContent(ipos+first)-baseline); + signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline); his->SetBinContent(ipos+5,signal); } his->Write(hname); @@ -233,7 +535,7 @@ void AliTPCCalibTCF::ProcessRawEvent(AliTPCRawStream *rawStream, const char *nam his->AddBinContent(1,1); // pulse counter for each pad for (Int_t ipos=0; iposGetBinContent(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); @@ -243,6 +545,8 @@ void AliTPCCalibTCF::ProcessRawEvent(AliTPCRawStream *rawStream, const char *nam tempHis->Reset(); } prevTime = time; + prevSec = sector; + prevRow = row; prevPad = pad; } @@ -274,7 +578,7 @@ void AliTPCCalibTCF::MergeHistoPerSector(const char *nameFileIn) { TIter next( fileIn.GetListOfKeys() ); char nameFileOut[100]; - sprintf(nameFileOut,"Sec-%s",nameFileIn); + snprintf(nameFileOut,100,"Sec-%s",nameFileIn); TFile fileOut(nameFileOut,"RECREATE"); fileOut.cd(); @@ -285,15 +589,18 @@ void AliTPCCalibTCF::MergeHistoPerSector(const char *nameFileIn) { 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) @@ -312,26 +619,24 @@ void AliTPCCalibTCF::MergeHistoPerSector(const char *nameFileIn) { for (Int_t ipos=0; iposAddBinContent(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 @@ -348,7 +653,7 @@ void AliTPCCalibTCF::AnalyzeRootFile(const char *nameFileIn, Int_t minNumPulse) TIter next( fileIn.GetListOfKeys() ); char nameFileOut[100]; - sprintf(nameFileOut,"TCFparam-%s",nameFileIn); + snprintf(nameFileOut,100,"TCF-%s",nameFileIn); TFile fileOut(nameFileOut,"RECREATE"); fileOut.cd(); @@ -359,13 +664,17 @@ void AliTPCCalibTCF::AnalyzeRootFile(const char *nameFileIn, Int_t minNumPulse) 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 + ++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++){ @@ -382,8 +691,10 @@ void AliTPCCalibTCF::AnalyzeRootFile(const char *nameFileIn, Int_t minNumPulse) } coefP->~Double_t(); coefZ->~Double_t(); + } else { + printf("Skip histogram %d out of %d | not enough accumulated pulses\n",iHist,nHist); } - + } fileIn.Close(); @@ -394,7 +705,7 @@ void AliTPCCalibTCF::AnalyzeRootFile(const char *nameFileIn, Int_t minNumPulse) //____________________________________________________________________________ -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 @@ -402,7 +713,7 @@ Int_t AliTPCCalibTCF::AnalyzePulse(TH1F *hisIn, Double_t *coefZ, Double_t *coefP // 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); @@ -430,8 +741,15 @@ Int_t AliTPCCalibTCF::AnalyzePulse(TH1F *hisIn, Double_t *coefZ, Double_t *coefP 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); @@ -455,7 +773,7 @@ Int_t AliTPCCalibTCF::AnalyzePulse(TH1F *hisIn, Double_t *coefZ, Double_t *coefP //____________________________________________________________________________ -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 @@ -477,7 +795,7 @@ void AliTPCCalibTCF::TestTCFonRootFile(const char *nameFileIn, const char *nameF } 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"); @@ -493,12 +811,16 @@ void AliTPCCalibTCF::TestTCFonRootFile(const char *nameFileIn, const char *nameF 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); @@ -525,7 +847,7 @@ void AliTPCCalibTCF::TestTCFonRootFile(const char *nameFileIn, const char *nameF //_____________________________________________________________________________ -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) @@ -540,9 +862,27 @@ void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameF // 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++){ @@ -550,44 +890,67 @@ void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameF 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); } @@ -598,7 +961,7 @@ void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameF 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.?) @@ -607,8 +970,10 @@ void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameF 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(); @@ -617,19 +982,41 @@ void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameF 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 && maxlowLim && max= 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; iposGetBinContent(ipos+first)-baseline); + signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline); his->SetBinContent(ipos+5,signal); } @@ -639,7 +1026,7 @@ void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameF // (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(); @@ -649,22 +1036,27 @@ void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameF 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(); @@ -673,35 +1065,130 @@ void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameF } +//____________________________________________________________________________ +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; rowiGetNRows(sec); rowi++) { + for (UInt_t padi=0; padiGetNPads(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]; @@ -709,38 +1196,37 @@ TNtuple *AliTPCCalibTCF::PlotOccupSummary(const char *nameFile, Int_t nPulseMin) 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(); } //____________________________________________________________________________ @@ -764,70 +1250,30 @@ void AliTPCCalibTCF::PlotQualitySummary(const char *nameFileQuality, const char 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; iGetEntries(); 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) { @@ -836,18 +1282,18 @@ Int_t AliTPCCalibTCF::FitPulse(TNtuple *dataTuple, Double_t *coefZ, Double_t *co // // 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 @@ -857,31 +1303,31 @@ Int_t AliTPCCalibTCF::FitPulse(TNtuple *dataTuple, Double_t *coefZ, Double_t *co 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 @@ -890,7 +1336,7 @@ Int_t AliTPCCalibTCF::FitPulse(TNtuple *dataTuple, Double_t *coefZ, Double_t *co 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; } @@ -906,7 +1352,7 @@ Int_t AliTPCCalibTCF::FitPulse(TNtuple *dataTuple, Double_t *coefZ, Double_t *co fitParam->~Double_t(); valuePZ->~Double_t(); - gMinuit->~TMinuit(); + minuitFit->~TMinuit(); return 1; @@ -916,7 +1362,7 @@ Int_t AliTPCCalibTCF::FitPulse(TNtuple *dataTuple, Double_t *coefZ, Double_t *co //____________________________________________________________________________ -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 @@ -980,10 +1426,10 @@ Double_t* AliTPCCalibTCF::ExtractPZValues(Double_t *param) { // 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]; @@ -991,7 +1437,7 @@ Double_t* AliTPCCalibTCF::ExtractPZValues(Double_t *param) { 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]; @@ -999,7 +1445,7 @@ Double_t* AliTPCCalibTCF::ExtractPZValues(Double_t *param) { 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]; @@ -1037,18 +1483,23 @@ Double_t* AliTPCCalibTCF::ExtractPZValues(Double_t *param) { //____________________________________________________________________________ -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; iposGetEntry(ipos); Float_t *p = dataTuple->GetArgs(); @@ -1088,20 +1539,29 @@ void AliTPCCalibTCF::Equalization(TNtuple *dataTuple, Double_t *coefZ, Double_t 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' @@ -1122,10 +1582,10 @@ Int_t AliTPCCalibTCF::FindCorTCFparam(TH1F *hisIn, const char *nameFileTCF, Doub 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! @@ -1137,9 +1597,9 @@ Int_t AliTPCCalibTCF::FindCorTCFparam(TH1F *hisIn, const char *nameFileTCF, Doub 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]; @@ -1150,8 +1610,8 @@ Int_t AliTPCCalibTCF::FindCorTCFparam(TH1F *hisIn, const char *nameFileTCF, Doub } 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; @@ -1233,7 +1693,7 @@ Double_t *AliTPCCalibTCF::GetQualityOfTCF(TH1F *hisIn, Double_t *coefZ, Double_t startOfPulse = 1; posOfStart = i; } - if( (sig < threshADC) && (startOfPulse == 1) ){ + if( (sig <= threshADC) && (startOfPulse == 1) ){ widthFound = 1; width = i - posOfStart + 1; } @@ -1244,7 +1704,7 @@ Double_t *AliTPCCalibTCF::GetQualityOfTCF(TH1F *hisIn, Double_t *coefZ, Double_t startOfPulseTCF = 1; posOfStartTCF = i; } - if( (sigTCF < threshADC) && (startOfPulseTCF == 1) ){ + if( (sigTCF <= threshADC) && (startOfPulseTCF == 1) ){ widthFoundTCF = 1; widthTCF = i -posOfStartTCF + 1; } @@ -1262,7 +1722,7 @@ Double_t *AliTPCCalibTCF::GetQualityOfTCF(TH1F *hisIn, Double_t *coefZ, Double_t } // 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; } } @@ -1294,7 +1754,7 @@ Double_t *AliTPCCalibTCF::GetQualityOfTCF(TH1F *hisIn, Double_t *coefZ, Double_t 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(); @@ -1307,7 +1767,7 @@ Double_t *AliTPCCalibTCF::GetQualityOfTCF(TH1F *hisIn, Double_t *coefZ, Double_t 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); } @@ -1326,7 +1786,7 @@ Double_t *AliTPCCalibTCF::GetQualityOfTCF(TH1F *hisIn, Double_t *coefZ, Double_t //____________________________________________________________________________ -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 @@ -1352,17 +1812,17 @@ TNtuple *AliTPCCalibTCF::ApplyTCFilter(TH1F *hisIn, Double_t *coefZ, Double_t *c } // 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]*(TMath::Power(2,16)-1)); - valZ[i] = (Int_t)(coefZ[i]*(TMath::Power(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; @@ -1393,7 +1853,7 @@ TNtuple *AliTPCCalibTCF::ApplyTCFilter(TH1F *hisIn, Double_t *coefZ, Double_t *c 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"); @@ -1403,21 +1863,16 @@ TNtuple *AliTPCCalibTCF::ApplyTCFilter(TH1F *hisIn, Double_t *coefZ, Double_t *c // 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() { // @@ -1430,5 +1885,394 @@ 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; iiAddBinContent(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; rowGetNpadrows(roc); row++) { + for (pad = 0; padGetNpads(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; + +}