// Origin: Marian Ivanov
//-------------------------------------------------------
-#include "AliTPCReconstructor.h"
-#include "AliTPCclustererMI.h"
-#include "AliTPCclusterMI.h"
-#include "AliTPCclusterInfo.h"
-#include <TObjArray.h>
+#include "Riostream.h"
+#include <TF1.h>
#include <TFile.h>
-#include "TGraph.h"
-#include "TF1.h"
-#include "TRandom.h"
-#include "AliMathBase.h"
+#include <TGraph.h>
+#include <TH1F.h>
+#include <TObjArray.h>
+#include <TRandom.h>
+#include <TTree.h>
+#include <TTreeStream.h>
+#include <TVirtualFFT.h>
-#include "AliTPCClustersArray.h"
-#include "AliTPCClustersRow.h"
#include "AliDigits.h"
-#include "AliSimDigits.h"
-#include "AliTPCParam.h"
-#include "AliTPCRecoParam.h"
-#include "AliRawReader.h"
-#include "AliTPCRawStream.h"
+#include "AliLoader.h"
+#include "AliLog.h"
+#include "AliMathBase.h"
#include "AliRawEventHeaderBase.h"
+#include "AliRawReader.h"
#include "AliRunLoader.h"
-#include "AliLoader.h"
-#include "Riostream.h"
-#include <TTree.h>
-#include "AliTPCcalibDB.h"
+#include "AliSimDigits.h"
#include "AliTPCCalPad.h"
#include "AliTPCCalROC.h"
-#include "TTreeStream.h"
-#include "AliLog.h"
-#include "TVirtualFFT.h"
+#include "AliTPCClustersArray.h"
+#include "AliTPCClustersRow.h"
+#include "AliTPCParam.h"
+#include "AliTPCRawStream.h"
+#include "AliTPCRecoParam.h"
+#include "AliTPCReconstructor.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCclusterInfo.h"
+#include "AliTPCclusterMI.h"
+#include "AliTPCclustererMI.h"
ClassImp(AliTPCclustererMI)
AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam):
fBins(0),
- fResBins(0),
+ fSigBins(0),
+ fNSigBins(0),
fLoop(0),
fMaxBin(0),
fMaxTime(0),
fAmplitudeHisto(0),
fDebugStreamer(0),
fRecoParam(0),
+ fBDumpSignal(kFALSE),
fFFTr2c(0)
{
//
AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI ¶m)
:TObject(param),
fBins(0),
- fResBins(0),
+ fSigBins(0),
+ fNSigBins(0),
fLoop(0),
fMaxBin(0),
fMaxTime(0),
fNcluster(0),
fAmplitudeHisto(0),
fDebugStreamer(0),
- fRecoParam(0)
+ fRecoParam(0),
+ fBDumpSignal(kFALSE),
+ fFFTr2c(0)
{
//
// dummy
//sigma z2 = in digits - angle estimated supposing vertex constraint
Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
- Float_t angular = fPadLength*(fParam->GetZLength()-z)/(fRx*fZWidth);
+ Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
angular*=angular;
angular/=12.;
Float_t sres = fParam->GetZSigma()/fZWidth;
// set pointers to data
//Int_t dummy[5] ={0,0,0,0,0};
Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2
- Float_t * resmatrix[5];
for (Int_t di=-2;di<=2;di++){
matrix[di+2] = &bins[k+di*max];
- resmatrix[di+2] = &fResBins[k+di*max];
}
//build matrix with virtual charge
Float_t sigmay2= GetSigmaY2(j0);
c.SetSigmaY2(mi2);
c.SetSigmaZ2(mj2);
AddCluster(c,(Float_t*)vmatrix,k);
- //remove cluster data from data
- for (Int_t di=-2;di<=2;di++)
- for (Int_t dj=-2;dj<=2;dj++){
- resmatrix[di+2][dj] -= vmatrix[di+2][dj+2];
- if (resmatrix[di+2][dj]<0) resmatrix[di+2][dj]=0;
- }
- resmatrix[2][0] =0;
return;
}
}
c.SetZ(fZWidth*(meanj-3));
c.SetZ(c.GetZ() - 3.*fParam->GetZSigma() + fParam->GetNTBinsL1()*fParam->GetZWidth()); // PASA delay + L1 delay
- c.SetZ(fSign*(fParam->GetZLength() - c.GetZ()));
+ c.SetZ(fSign*(fParam->GetZLength(fSector) - c.GetZ()));
c.SetX(fRx);
c.SetDetector(fSector);
c.SetRow(fRow);
TClonesArray * arr = fRowCl->GetArray();
AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
- if (matrix) {
+ if (fRecoParam->DumpSignal() &&matrix ) {
Int_t nbins=0;
Float_t *graph =0;
- if (fRecoParam->GetCalcPedestal()){
+ if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
nbins = fMaxTime;
graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
}
AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
cl->SetInfo(info);
}
+ if (!fRecoParam->DumpSignal()) {
+ cl->SetInfo(0);
+ }
fNcluster++;
}
fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
fBins =new Float_t[fMaxBin];
- fResBins =new Float_t[fMaxBin]; //fBins with residuals after 1 finder loop
+ fSigBins =new Int_t[fMaxBin];
+ fNSigBins = 0;
memset(fBins,0,sizeof(Float_t)*fMaxBin);
- memset(fResBins,0,sizeof(Float_t)*fMaxBin);
if (digarr.First()) //MI change
do {
if (dig<=fParam->GetZeroSup()) continue;
Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
- fBins[i*fMaxTime+j]=dig/gain;
+ Int_t bin = i*fMaxTime+j;
+ fBins[bin]=dig/gain;
+ fSigBins[fNSigBins++]=bin;
} while (digarr.Next());
digarr.ExpandTrackBuffer();
fOutput->Fill();
delete clrow;
nclusters+=fNcluster;
- delete[] fBins;
- delete[] fResBins;
+ delete[] fBins;
+ delete[] fSigBins;
}
Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
//alocate memory for sector - maximal case
//
Float_t** allBins = NULL;
- Float_t** allBinsRes = NULL;
+ Int_t** allSigBins = NULL;
+ Int_t* allNSigBins = NULL;
Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
allBins = new Float_t*[nRowsMax];
- allBinsRes = new Float_t*[nRowsMax];
+ allSigBins = new Int_t*[nRowsMax];
+ allNSigBins = new Int_t[nRowsMax];
for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
//
Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
allBins[iRow] = new Float_t[maxBin];
- allBinsRes[iRow] = new Float_t[maxBin];
memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
+ allSigBins[iRow] = new Int_t[maxBin];
+ allNSigBins[iRow]=0;
}
//
// Loop over sectors
Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
+ allNSigBins[iRow] = 0;
}
// Loas the raw data for corresponding DDLs
Float_t signal = input.GetSignal();
if (!calcPedestal && signal <= zeroSup) continue;
if (!calcPedestal) {
- allBins[iRow][iPad*fMaxTime+iTimeBin] = signal/gain;
+ Int_t bin = iPad*fMaxTime+iTimeBin;
+ allBins[iRow][bin] = signal/gain;
+ allSigBins[iRow][allNSigBins[iRow]++] = bin;
}else{
allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
}
//
for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
- allBins[iRow][iPad*fMaxTime+iTimeBin] -= pedestalEvent;
+ Int_t bin = iPad*fMaxTime+iTimeBin;
+ allBins[iRow][bin] -= pedestalEvent;
if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
- allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
+ allBins[iRow][bin] = 0;
if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
- allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
+ allBins[iRow][bin] = 0;
if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
- allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
- if (allBins[iRow][iPad*fMaxTime+iTimeBin] < 3.0*rmsEvent) // 3 sigma cut on RMS
- allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
+ allBins[iRow][bin] = 0;
+ if (allBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
+ allBins[iRow][bin] = 0;
+ if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
}
}
}
fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
fBins = allBins[fRow];
- fResBins = allBinsRes[fRow];
+ fSigBins = allSigBins[fRow];
+ fNSigBins = allNSigBins[fRow];
FindClusters(noiseROC);
for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
delete [] allBins[iRow];
- delete [] allBinsRes[iRow];
+ delete [] allSigBins[iRow];
}
delete [] allBins;
- delete [] allBinsRes;
+ delete [] allSigBins;
+ delete [] allNSigBins;
Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
//
// add virtual charge at the edge
//
- if (0) for (Int_t i=0; i<fMaxTime; i++){
- Float_t amp1 = fBins[i+3*fMaxTime];
- Float_t amp0 =0;
- if (amp1>0){
- Float_t amp2 = fBins[i+4*fMaxTime];
- if (amp2==0) amp2=0.5;
- Float_t sigma2 = GetSigmaY2(i);
- amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);
- if (gDebug>4) printf("\n%f\n",amp0);
- }
- fBins[i+2*fMaxTime] = amp0;
- amp0 = 0;
- amp1 = fBins[(fMaxPad+2)*fMaxTime+i];
- if (amp1>0){
- Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime];
- if (amp2==0) amp2=0.5;
- Float_t sigma2 = GetSigmaY2(i);
- amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);
- if (gDebug>4) printf("\n%f\n",amp0);
- }
- fBins[(fMaxPad+3)*fMaxTime+i] = amp0;
- }
- memcpy(fResBins,fBins, fMaxBin*sizeof(Float_t));
- //
- //
+ Double_t kMaxDumpSize = 500000;
+ if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
//
fNcluster=0;
fLoop=1;
- Float_t *b=&fBins[-1]+2*fMaxTime;
- Int_t crtime = Int_t((fParam->GetZLength()-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
+ Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
- for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
- b++;
- if (i%fMaxTime<crtime) {
- Int_t delta = -(i%fMaxTime)+crtime;
- b+=delta;
- i+=delta;
- continue;
- }
+ for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
+ Int_t i = fSigBins[iSig];
+ if (i%fMaxTime<=crtime) continue;
+ Float_t *b = &fBins[i];
//absolute custs
if (b[0]<minMaxCutAbs) continue; //threshold for maxima
//
// ESTIMATE pedestal and the noise
//
const Int_t kPedMax = 100;
+ Double_t kMaxDebugSize = 5000000.;
Float_t max = 0;
Float_t maxPos = 0;
Int_t median = -1;
}
// truncated mean
//
- Double_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
- Double_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
- Double_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
+ Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
+ Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
+ Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
//
for (Int_t idelta=1; idelta<10; idelta++){
if (median-idelta<=0) continue;
//
// Digital noise
//
- if (max-median>30.*TMath::Max(1.,rms06)){
- //
- //
- TGraph * graph =new TGraph(nchannels, dtime, dsignal);
- //
- //
- // jumps left - right
- Int_t njumps0=0;
- Double_t deltaT0[2000];
- Double_t deltaA0[2000];
- Int_t lastJump0 = fRecoParam->GetFirstBin();
- Int_t njumps1=0;
- Double_t deltaT1[2000];
- Double_t deltaA1[2000];
- Int_t lastJump1 = fRecoParam->GetFirstBin();
- Int_t njumps2=0;
- Double_t deltaT2[2000];
- Double_t deltaA2[2000];
- Int_t lastJump2 = fRecoParam->GetFirstBin();
-
- for (Int_t itime=fRecoParam->GetFirstBin()+1; itime<fRecoParam->GetLastBin()-1; itime++){
- if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,rms06) &&
- TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,rms06) &&
- (dsignal[itime-1]-median<5.*rms06) &&
- (dsignal[itime+1]-median<5.*rms06)
- ){
- deltaA0[njumps0] = dsignal[itime]-dsignal[itime-1];
- deltaT0[njumps0] = itime-lastJump0;
- lastJump0 = itime;
- njumps0++;
- }
- if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,rms06) &&
- (dsignal[itime-1]-median<5.*rms06)
- ) {
- deltaA1[njumps1] = dsignal[itime]-dsignal[itime-1];
- deltaT1[njumps1] = itime-lastJump1;
- lastJump1 = itime;
- njumps1++;
- }
- if (TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,rms06) &&
- (dsignal[itime+1]-median<5.*rms06)
- ) {
- deltaA2[njumps2] = dsignal[itime]-dsignal[itime+1];
- deltaT2[njumps2] = itime-lastJump2;
- lastJump2 = itime;
- njumps2++;
- }
- }
- //
- if (njumps0>0 || njumps1>0 || njumps2>0){
- TGraph *graphDN0 = new TGraph(njumps0, deltaT0, deltaA0);
- TGraph *graphDN1 = new TGraph(njumps1, deltaT1, deltaA1);
- TGraph *graphDN2 = new TGraph(njumps2, deltaT2, deltaA2);
- (*fDebugStreamer)<<"SignalDN"<< //digital - noise pads - or random sample of pads
- "TimeStamp="<<fTimeStamp<<
- "EventType="<<fEventType<<
- "Sector="<<uid[0]<<
- "Row="<<uid[1]<<
- "Pad="<<uid[2]<<
- "Graph="<<graph<<
- "Max="<<max<<
- "MaxPos="<<maxPos<<
- "Graph.="<<graph<<
- "P0GraphDN0.="<<graphDN0<<
- "P1GraphDN1.="<<graphDN1<<
- "P2GraphDN2.="<<graphDN2<<
- //
- "Median="<<median<<
- "Mean="<<mean<<
- "RMS="<<rms<<
- "Mean06="<<mean06<<
- "RMS06="<<rms06<<
- "Mean09="<<mean09<<
- "RMS09="<<rms09<<
- "\n";
- delete graphDN0;
- delete graphDN1;
- delete graphDN2;
- }
- delete graph;
- }
+ // if (max-median>30.*TMath::Max(1.,Double_t(rms06)) && (((*fDebugStreamer)<<"SignalDN").GetSize()<kMaxDebugSize)){
+// //
+// //
+// TGraph * graph =new TGraph(nchannels, dtime, dsignal);
+// //
+// //
+// // jumps left - right
+// Int_t njumps0=0;
+// Double_t deltaT0[2000];
+// Double_t deltaA0[2000];
+// Int_t lastJump0 = fRecoParam->GetFirstBin();
+// Int_t njumps1=0;
+// Double_t deltaT1[2000];
+// Double_t deltaA1[2000];
+// Int_t lastJump1 = fRecoParam->GetFirstBin();
+// Int_t njumps2=0;
+// Double_t deltaT2[2000];
+// Double_t deltaA2[2000];
+// Int_t lastJump2 = fRecoParam->GetFirstBin();
+
+// for (Int_t itime=fRecoParam->GetFirstBin()+1; itime<fRecoParam->GetLastBin()-1; itime++){
+// if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06)) &&
+// TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06)) &&
+// (dsignal[itime-1]-median<5.*rms06) &&
+// (dsignal[itime+1]-median<5.*rms06)
+// ){
+// deltaA0[njumps0] = dsignal[itime]-dsignal[itime-1];
+// deltaT0[njumps0] = itime-lastJump0;
+// lastJump0 = itime;
+// njumps0++;
+// }
+// if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06)) &&
+// (dsignal[itime-1]-median<5.*rms06)
+// ) {
+// deltaA1[njumps1] = dsignal[itime]-dsignal[itime-1];
+// deltaT1[njumps1] = itime-lastJump1;
+// lastJump1 = itime;
+// njumps1++;
+// }
+// if (TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06)) &&
+// (dsignal[itime+1]-median<5.*rms06)
+// ) {
+// deltaA2[njumps2] = dsignal[itime]-dsignal[itime+1];
+// deltaT2[njumps2] = itime-lastJump2;
+// lastJump2 = itime;
+// njumps2++;
+// }
+// }
+// //
+// if (njumps0>0 || njumps1>0 || njumps2>0){
+// TGraph *graphDN0 = new TGraph(njumps0, deltaT0, deltaA0);
+// TGraph *graphDN1 = new TGraph(njumps1, deltaT1, deltaA1);
+// TGraph *graphDN2 = new TGraph(njumps2, deltaT2, deltaA2);
+// (*fDebugStreamer)<<"SignalDN"<< //digital - noise pads - or random sample of pads
+// "TimeStamp="<<fTimeStamp<<
+// "EventType="<<fEventType<<
+// "Sector="<<uid[0]<<
+// "Row="<<uid[1]<<
+// "Pad="<<uid[2]<<
+// "Graph="<<graph<<
+// "Max="<<max<<
+// "MaxPos="<<maxPos<<
+// "Graph.="<<graph<<
+// "P0GraphDN0.="<<graphDN0<<
+// "P1GraphDN1.="<<graphDN1<<
+// "P2GraphDN2.="<<graphDN2<<
+// //
+// "Median="<<median<<
+// "Mean="<<mean<<
+// "RMS="<<rms<<
+// "Mean06="<<mean06<<
+// "RMS06="<<rms06<<
+// "Mean09="<<mean09<<
+// "RMS09="<<rms09<<
+// "\n";
+// delete graphDN0;
+// delete graphDN1;
+// delete graphDN2;
+// }
+// delete graph;
+// }
//
// NOISE STUDY Fourier transform
//
TGraph * graph;
Bool_t random = (gRandom->Rndm()<0.0003);
- if (max-median>kMin || rms06>1.*fParam->GetZeroSup() || random){
+ if (((*fDebugStreamer)<<"SignalN").GetSize()<kMaxDebugSize)
+ if (max-median>kMin || rms06>1.*fParam->GetZeroSup() || random){
graph =new TGraph(nchannels, dtime, dsignal);
if (rms06>1.*fParam->GetZeroSup() || random){
//Double_t *input, Double_t threshold, Bool_t locMax, Double_t *freq, Double_t *re, Double_t *im, Double_t *mag, Double_t *phi);
//
// Central Electrode signal analysis
//
- Double_t ceQmax =0, ceQsum=0, ceTime=0;
- Double_t cemean = mean06, cerms=rms06 ;
+ Float_t ceQmax =0, ceQsum=0, ceTime=0;
+ Float_t cemean = mean06, cerms=rms06 ;
Int_t cemaxpos= 0;
- Double_t ceThreshold=5.*cerms;
- Double_t ceSumThreshold=8.*cerms;
+ Float_t ceThreshold=5.*cerms;
+ Float_t ceSumThreshold=8.*cerms;
const Int_t kCemin=5; // range for the analysis of the ce signal +- channels from the peak
const Int_t kCemax=5;
for (Int_t i=nchannels-2; i>nchannels/2; i--){
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff