Debug information attached to AliTPCclusterMI in debug mode
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 22 Jan 2007 14:35:23 +0000 (14:35 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 22 Jan 2007 14:35:23 +0000 (14:35 +0000)
TPC/AliTPCclusterInfo.cxx [new file with mode: 0644]
TPC/AliTPCclusterInfo.h [new file with mode: 0644]
TPC/AliTPCclusterMI.cxx
TPC/AliTPCclusterMI.h
TPC/AliTPCclustererMI.cxx
TPC/AliTPCclustererMI.h

diff --git a/TPC/AliTPCclusterInfo.cxx b/TPC/AliTPCclusterInfo.cxx
new file mode 100644 (file)
index 0000000..7e5b04d
--- /dev/null
@@ -0,0 +1,141 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-------------------------------------------------------
+//          Implementation of the TPC cluster debug information
+//
+//   Origin: Marian Ivanov   Marian.Ivanov@cern.ch
+//   
+//   Additional cluster information to monitor clustering performance
+//   and to extract a features of detector response
+//   Information attached to the default cluster
+//   ONLY in DEBUG MODE
+//    
+//-------------------------------------------------------
+
+/* $Id$ */
+
+#include "AliTPCclusterInfo.h"
+#include "AliLog.h"
+
+ClassImp(AliTPCclusterInfo);
+
+
+AliTPCclusterInfo::AliTPCclusterInfo():
+  fNPads(0),
+  fNTimeBins(0),
+  fNBins(0),
+  fGraph(0)
+{
+  //
+  // default constructor
+  //  
+  for (Int_t i=0; i<25;i++){
+    fMatrix[i] = i; 
+  }
+}
+
+AliTPCclusterInfo::AliTPCclusterInfo(const  AliTPCclusterInfo & info):
+  TObject(info),
+  fNPads(info.fNPads),
+  fNTimeBins(info.fNTimeBins),
+  fNBins(info.fNBins),
+  fGraph(0)    
+{
+  //
+  // copy constructor
+  //
+  //  AliInfo("Copy constructor\n");
+  for (Int_t i=0; i<25;i++){
+    fMatrix[i] = info.fMatrix[i]; 
+  }
+  if (info.fGraph) fGraph = new Float_t[fNBins];
+  for (Int_t i=0;i<fNBins; i++){
+    fGraph[i] = info.fGraph[i];
+  }
+  
+}
+
+
+AliTPCclusterInfo::AliTPCclusterInfo(Bool_t extend):
+  fNPads(0),
+  fNTimeBins(0),
+  fNBins(0),
+  fGraph(0)
+{
+  //
+  // allocate dummy graph - neccessary for IO part to use automatic branching
+  //  
+  for (Int_t i=0; i<25;i++){
+    fMatrix[i] = i; 
+  }
+  if (extend){
+    fNBins = 1;
+    fGraph  = new Float_t[1];
+    fGraph[0]=-1;
+  }
+}
+
+AliTPCclusterInfo::AliTPCclusterInfo(Float_t *matrix, Int_t nbins, Float_t* graph){
+  //
+  // constructor of the info
+  //
+  for (Int_t i=0;i<25;i++){
+    fMatrix[i]=matrix[i];
+  }
+  fNPads=0;
+  fNTimeBins=0;
+  Int_t center = 5+5+2;
+  for (Int_t i=-2; i<=2;i++) if (matrix[center+i]>0) fNTimeBins++;
+  for (Int_t i=-2; i<=2;i++) if (matrix[center+i*5]>0) fNPads++;
+  fNBins = nbins;
+  fGraph = 0;
+  if (fNBins>0) {
+    fGraph = new Float_t[fNBins];
+    for (Int_t i=0;i<fNBins; i++){
+      fGraph[i] = graph[i];
+    }
+  }
+}
+
+UChar_t AliTPCclusterInfo::GetNPads(Float_t threshold) const { 
+  //
+  //
+  Int_t nPads=0;
+  Int_t center = 5+5+2;
+  for (Int_t i=-2; i<=2;i++) if (fMatrix[center+i*5]>threshold) nPads++;
+  return nPads;
+}
+
+UChar_t AliTPCclusterInfo::GetNTimeBins(Float_t threshold) const { 
+  //
+  //
+  //
+  Int_t nTimeBins=0;
+  Int_t center = 5+5+2;
+  for (Int_t i=-2; i<=2;i++) if (fMatrix[center+i]>0) nTimeBins++;
+  return nTimeBins;
+}
+
+
+
+
+AliTPCclusterInfo::~AliTPCclusterInfo(){
+  //
+  // destructor
+  //
+  if (fGraph)  delete [] fGraph;
+}
+
diff --git a/TPC/AliTPCclusterInfo.h b/TPC/AliTPCclusterInfo.h
new file mode 100644 (file)
index 0000000..e3df089
--- /dev/null
@@ -0,0 +1,44 @@
+#ifndef ALITPCCLUSTERINFO_H
+#define ALITPCCLUSTERINFO_H
+
+//-------------------------------------------------------
+//                    TPC Cluster Class
+//   Information for debugging puposes
+//   Origin: Marian Ivanov
+//-------------------------------------------------------
+
+/* $Id$ */
+
+
+#include "TObject.h"
+//_____________________________________________________________________________
+class AliTPCclusterInfo : public TObject {
+public:
+  AliTPCclusterInfo();
+  AliTPCclusterInfo(Bool_t withGraph);
+  AliTPCclusterInfo(Float_t *matrix, Int_t nbins, Float_t* graph);
+  AliTPCclusterInfo(const  AliTPCclusterInfo & info);
+  virtual ~AliTPCclusterInfo();
+  UChar_t GetNPads() const { return fNPads;}
+  UChar_t GetNTimeBins() const { return fNTimeBins;}
+  UChar_t GetNPads(Float_t threshold) const;
+  UChar_t GetNTimeBins(Float_t threshold) const;
+  Float_t* GetMatrix(){ return fMatrix;}
+  void    SetGraph(Float_t * graph, Int_t nbins){ fGraph = graph; fNBins=nbins;}
+  void    SetNPadsTimes(UChar_t npads, UChar_t ntimes){ fNPads = npads; fNTimeBins = ntimes;}
+ protected:
+  Float_t  fMatrix[25];               // matrix of amplitude arround center pad - time 
+  UChar_t  fNPads;                // number of pads      in cluster
+  UChar_t  fNTimeBins;            // number of time bins in cluster
+  Int_t  fNBins;                // number of bins in graph
+  Float_t *fGraph;               //[fNBins] signal time dependence graph
+private:
+  ClassDef(AliTPCclusterInfo,1)  // Time Projection Chamber clusters Inofrmation 
+};
+
+
+
+
+#endif
+
+
index 0b8194b..efd550d 100644 (file)
 /* $Id$ */
 
 #include "AliTPCclusterMI.h"
+#include "AliTPCclusterInfo.h"
+#include "AliLog.h"
 
 ClassImp(AliTPCclusterMI)
 
 
-AliTPCclusterMI::AliTPCclusterMI():
+AliTPCclusterMI::AliTPCclusterMI(Bool_t withInfo):
   AliCluster(),
   fX(0),
   fQ(0),
@@ -39,13 +41,61 @@ AliTPCclusterMI::AliTPCclusterMI():
   fMax(0),
   fUsed(0),
   fDetector(0),
-  fRow(0)
+  fRow(0),
+  fTimeBin(0),
+  fPad(0),
+  fInfo(0)
 {
   //
   // default constructor
   //
+  if (withInfo) fInfo = new AliTPCclusterInfo;
 }
 
+AliTPCclusterMI::AliTPCclusterMI(const AliTPCclusterMI & cluster):
+  AliCluster(cluster),
+  fX(cluster.fX),
+  fQ(cluster.fQ),
+  fType(cluster.fType),
+  fMax(cluster.fMax),
+  fUsed(cluster.fUsed),
+  fDetector(cluster.fDetector),
+  fRow(cluster.fRow),
+  fTimeBin(cluster.fTimeBin),
+  fPad(cluster.fPad),
+  fInfo(0)
+{
+  //
+  // copy constructor
+  // 
+  //  AliInfo("Copy constructor\n");
+  if (cluster.fInfo) fInfo = new AliTPCclusterInfo(*(cluster.fInfo));
+}
+
+AliTPCclusterMI & AliTPCclusterMI::operator = (const AliTPCclusterMI & cluster)
+{
+  //
+  // assignment operator
+  // 
+  //  AliInfo("Asignment operator\n");
+
+  (AliCluster&)(*this) = (AliCluster&)cluster;
+  fX    = cluster.fX;
+  fQ    = cluster.fQ;
+  fType = cluster.fType;
+  fMax  = cluster.fMax;
+  fUsed = cluster.fUsed;
+  fDetector = cluster.fDetector;
+  fRow  = cluster.fRow;
+  fTimeBin = cluster.fTimeBin;
+  fPad     = cluster.fPad;
+  fInfo = 0;
+  if (cluster.fInfo) fInfo = new AliTPCclusterInfo(*(cluster.fInfo));
+}
+
+
+
+
 AliTPCclusterMI::AliTPCclusterMI(Int_t *lab, Float_t *hit) : 
   AliCluster(lab,hit),
   fX(0),
@@ -54,14 +104,25 @@ AliTPCclusterMI::AliTPCclusterMI(Int_t *lab, Float_t *hit) :
   fMax(0),
   fUsed(0),
   fDetector(0),
-  fRow(0)    
+  fRow(0),
+  fInfo(0)
 {
   //
   // constructor
   //
   fQ = (UShort_t)hit[4];
+  fInfo = 0;
 }
 
+AliTPCclusterMI::~AliTPCclusterMI() {
+  //
+  // destructor
+  //
+  if (fInfo) delete fInfo;
+}
+
+
+
 Bool_t AliTPCclusterMI::IsSortable() const
 {
   //
index dc3700c..421f44c 100644 (file)
 
 #include "AliCluster.h"
 #include "TMath.h"
-
+class AliTPCclusterInfo;
 //_____________________________________________________________________________
 class AliTPCclusterMI : public AliCluster {
 public:
-  AliTPCclusterMI();
+  AliTPCclusterMI(Bool_t withInfo=kTRUE);
+  AliTPCclusterMI(const AliTPCclusterMI & cluster);
+  AliTPCclusterMI &operator = (const AliTPCclusterMI & cluster); //assignment operator
   AliTPCclusterMI(Int_t *lab, Float_t *hit);
-  virtual ~AliTPCclusterMI() {}
+  virtual ~AliTPCclusterMI();
   virtual Bool_t IsSortable() const; 
   virtual Int_t Compare(const TObject* obj) const;
   inline  void Use(Int_t inc=10);
@@ -28,6 +30,8 @@ public:
   virtual Int_t GetRow() const {return fRow;}
   virtual void SetDetector(Int_t detector){fDetector = (UChar_t)(detector%256);}
   virtual void SetRow(Int_t row){fRow = (UChar_t)(row%256);}  
+  virtual void SetTimeBin(Float_t timeBin){ fTimeBin= timeBin;}
+  virtual void SetPad(Float_t pad){ fPad = pad;}
   //
   void SetQ(Float_t q) {fQ=(UShort_t)q;}
   void SetType(Char_t type) {fType=type;}
@@ -36,7 +40,11 @@ public:
   Float_t GetQ() const {return TMath::Abs(fQ);}
   Float_t GetMax() const {return fMax;} 
   Char_t  GetType()const {return fType;}
+  Float_t GetTimeBin() const { return fTimeBin;}
+  Float_t GetPad() const { return fPad;}
+  AliTPCclusterInfo * GetInfo() const { return fInfo;}
+  void SetInfo(AliTPCclusterInfo * info) { fInfo = info;}
+
 private:
   Float_t   fX;        //X position of cluster
   Short_t   fQ ;       //Q of cluster (in ADC counts)  
@@ -45,7 +53,10 @@ private:
   Char_t    fUsed;     //counter of usage  
   UChar_t   fDetector; //detector  number
   UChar_t   fRow;      //row number number
-  ClassDef(AliTPCclusterMI,2)  // Time Projection Chamber clusters
+  Float_t   fTimeBin;  //time bin coordinate
+  Float_t   fPad;  //pad coordinate
+  AliTPCclusterInfo * fInfo;  // pointer to the cluster debug info
+  ClassDef(AliTPCclusterMI,3)  // Time Projection Chamber clusters
 };
 
 void AliTPCclusterMI::Use(Int_t inc) 
index 82d5ee2..b132ab1 100644 (file)
@@ -24,6 +24,7 @@
 #include "AliTPCReconstructor.h"
 #include "AliTPCclustererMI.h"
 #include "AliTPCclusterMI.h"
+#include "AliTPCclusterInfo.h"
 #include <TObjArray.h>
 #include <TFile.h>
 #include "TGraph.h"
@@ -49,7 +50,7 @@
 #include "AliTPCCalROC.h"
 #include "TTreeStream.h"
 #include "AliLog.h"
-
+#include "TVirtualFFT.h"
 
 ClassImp(AliTPCclustererMI)
 
@@ -82,7 +83,8 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar
   fNcluster(0),
   fAmplitudeHisto(0),
   fDebugStreamer(0),
-  fRecoParam(0)
+  fRecoParam(0),
+  fFFTr2c(0)
 {
   //
   // COSNTRUCTOR
@@ -102,6 +104,8 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar
   }
   fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
   fAmplitudeHisto = 0;
+  Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
+  fFFTr2c = TVirtualFFT::FFT(1, &nPoints, "R2C  K");
 }
 //______________________________________________________________
 AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
@@ -313,9 +317,11 @@ AliTPCclusterMI &c)
     c.SetQ(sumw);
     c.SetY(meani*fPadWidth); 
     c.SetZ(meanj*fZWidth); 
+    c.SetPad(meani);
+    c.SetTimeBin(meanj);
     c.SetSigmaY2(mi2);
     c.SetSigmaZ2(mj2);
-    AddCluster(c);
+    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++){
@@ -349,10 +355,12 @@ AliTPCclusterMI &c)
   c.SetQ(sumu);
   c.SetY(meani*fPadWidth); 
   c.SetZ(meanj*fZWidth); 
+  c.SetPad(meani);
+  c.SetTimeBin(meanj);
   c.SetSigmaY2(mi2);
   c.SetSigmaZ2(mj2);
   c.SetType(Char_t(overlap)+1);
-  AddCluster(c);
+  AddCluster(c,(Float_t*)vmatrix,k);
 
   //unfolding 2
   meani-=i0;
@@ -510,7 +518,7 @@ Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, F
   return max;
 }
 
-void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){
+void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){
   //
   // transform cluster to the global coordinata
   // add the cluster to the array
@@ -560,8 +568,17 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){
   if (fLoop==2) c.SetType(100);
 
   TClonesArray * arr = fRowCl->GetArray();
-  // AliTPCclusterMI * cl = 
-  new ((*arr)[fNcluster]) AliTPCclusterMI(c);
+  AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
+  if (matrix) {
+    Int_t nbins=0;
+    Float_t *graph =0;
+    if (fRecoParam->GetCalcPedestal()){
+      nbins = fMaxTime;
+      graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
+    }
+    AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
+    cl->SetInfo(info);
+  }
 
   fNcluster++;
 }
@@ -922,10 +939,13 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
       continue; 
     }
     //absolute custs
-    if (b[0]<minMaxCutAbs) continue;   //threshold form maxima  
-    if (b[-1]+b[1]<=0) continue;               // cut on isolated clusters
-    if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
-    if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue;   //threshold for up town TRF 
+    if (b[0]<minMaxCutAbs) continue;   //threshold for maxima  
+    //
+    if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue;  // cut on isolated clusters 
+    //    if (b[-1]+b[1]<=0) continue;               // cut on isolated clusters
+    //if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
+    //
+    if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue;   //threshold for up down  (TRF) 
     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue;   //threshold for left right (PRF)    
     if (!IsMaximum(*b,fMaxTime,b)) continue;
     //
@@ -935,7 +955,7 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
     if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue;   //threshold for up town TRF 
     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue;   //threshold for left right (PRF)    
   
-    AliTPCclusterMI c;
+    AliTPCclusterMI c(kFALSE);   // default cosntruction  without info
     Int_t dummy=0;
     MakeCluster(i, fMaxTime, fBins, dummy,c);
     
@@ -955,7 +975,6 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
   //
   // ESTIMATE pedestal and the noise
   // 
-  Int_t offset =100;
   const Int_t kPedMax = 100;
   Float_t  max    =  0;
   Float_t  maxPos =  0;
@@ -1095,19 +1114,115 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
     dsignal[i] = signal[i];
   }
   //
+  // 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;
+  }
+
+  //
+  // NOISE STUDY  Fourier transform
   //
   TGraph * graph;
-  Bool_t random = (gRandom->Rndm()<0.0001);
-  if (max-median>kMin || rms06>2.*fParam->GetZeroSup() || random){
+  Bool_t random = (gRandom->Rndm()<0.0003);
+  if (max-median>kMin || rms06>1.*fParam->GetZeroSup() || random){
     graph =new TGraph(nchannels, dtime, dsignal);
-    if (rms06>2.*fParam->GetZeroSup() || random)
+    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);
+      Float_t * input = &(dsignal[fRecoParam->GetFirstBin()]);
+      Float_t freq[2000], re[2000], im[2000], mag[2000], phi[2000];
+      Int_t npoints = TransformFFT(input, -1,kFALSE, freq, re, im, mag, phi);
+      TGraph *graphMag0 = new TGraph(npoints, freq, mag);
+      TGraph *graphPhi0 = new TGraph(npoints, freq, phi);
+      npoints = TransformFFT(input, 0.5,kTRUE, freq, re, im, mag, phi);
+      TGraph *graphMag1 = new TGraph(npoints, freq, mag);
+      TGraph *graphPhi1 = new TGraph(npoints, freq, phi);
+      
       (*fDebugStreamer)<<"SignalN"<<    //noise pads - or random sample of pads
        "TimeStamp="<<fTimeStamp<<
        "EventType="<<fEventType<<
        "Sector="<<uid[0]<<
        "Row="<<uid[1]<<
        "Pad="<<uid[2]<<
-       "Graph="<<graph<<
+       "Graph.="<<graph<<
        "Max="<<max<<
        "MaxPos="<<maxPos<<
        //
@@ -1118,7 +1233,21 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
        "RMS06="<<rms06<<
        "Mean09="<<mean09<<
        "RMS09="<<rms09<<
+       // FFT part
+       "Mag0.="<<graphMag0<<
+       "Mag1.="<<graphMag1<<
+       "Phi0.="<<graphPhi0<<
+       "Phi1.="<<graphPhi1<<
        "\n";
+      delete graphMag0;
+      delete graphMag1;
+      delete graphPhi0;
+      delete graphPhi1;
+    }
+    //
+    // Big signals dumping
+    //
+    
     if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin()) 
       (*fDebugStreamer)<<"SignalB"<<     // pads with signal
        "TimeStamp="<<fTimeStamp<<
@@ -1272,6 +1401,8 @@ void AliTPCclustererMI::DumpHistos(){
       Float_t rms  =  histo->GetRMS();
       Float_t gmean = histo->GetFunction("gaus")->GetParameter(1);
       Float_t gsigma = histo->GetFunction("gaus")->GetParameter(2);
+      Float_t gmeanErr = histo->GetFunction("gaus")->GetParError(1);
+      Float_t gsigmaErr = histo->GetFunction("gaus")->GetParError(2);
       Float_t max = histo->GetFunction("gaus")->GetParameter(0);
 
       // get pad number
@@ -1297,8 +1428,58 @@ void AliTPCclustererMI::DumpHistos(){
        "RMS="<<rms<<      
        "GMean="<<gmean<<
        "GSigma="<<gsigma<<
+       "GMeanErr="<<gmeanErr<<
+       "GSigmaErr="<<gsigmaErr<<
        "\n";
       if (array->UncheckedAt(ipad)) fDebugStreamer->StoreObject(array->UncheckedAt(ipad));
     }
   }
 }
+
+
+
+Int_t  AliTPCclustererMI::TransformFFT(Float_t *input, Float_t threshold, Bool_t locMax, Float_t *freq, Float_t *re, Float_t *im, Float_t *mag, Float_t *phi)
+{
+  //
+  // calculate fourrie transform 
+  // return only frequncies with mag over threshold
+  // if locMax is spectified only freque with local maxima over theshold is returned 
+
+  if (! fFFTr2c) return kFALSE;
+  if (!freq) return kFALSE;
+
+  Int_t current=0;
+  Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
+  Double_t *in = new Double_t[nPoints];
+  Double_t *rfft = new Double_t[nPoints];
+  Double_t *ifft = new Double_t[nPoints];
+  for (Int_t i=0; i<nPoints; i++){in[i]=input[i];}
+  fFFTr2c->SetPoints(in);
+  fFFTr2c->Transform();
+  fFFTr2c->GetPointsComplex(rfft, ifft);
+  for (Int_t i=3; i<nPoints/2-3; i++){
+    Float_t lmag =  TMath::Sqrt(rfft[i]*rfft[i]+ifft[i]*ifft[i])/nPoints;
+    if (lmag<threshold) continue;
+    if (locMax){
+      if ( TMath::Sqrt(rfft[i-1]*rfft[i-1]+ifft[i-1]*ifft[i-1])/nPoints>lmag) continue;
+      if ( TMath::Sqrt(rfft[i+1]*rfft[i+1]+ifft[i+1]*ifft[i+1])/nPoints>lmag) continue;
+      if ( TMath::Sqrt(rfft[i-2]*rfft[i-2]+ifft[i-2]*ifft[i-2])/nPoints>lmag) continue;
+      if ( TMath::Sqrt(rfft[i+2]*rfft[i+2]+ifft[i+2]*ifft[i+2])/nPoints>lmag) continue;
+      if ( TMath::Sqrt(rfft[i-3]*rfft[i-3]+ifft[i-3]*ifft[i-3])/nPoints>lmag) continue;
+      if ( TMath::Sqrt(rfft[i+3]*rfft[i+3]+ifft[i+3]*ifft[i+3])/nPoints>lmag) continue;
+    }
+    
+    freq[current] = Float_t(i)/Float_t(nPoints);
+    //
+    re[current] = rfft[i];
+    im[current] = ifft[i];
+    mag[current]=lmag;
+    phi[current]=TMath::ATan2(ifft[i],rfft[i]);
+    current++;
+  }
+  delete [] in;
+  delete [] rfft;
+  delete [] ifft;
+  return current;
+}
+
index 1fa57b1..3261140 100644 (file)
@@ -28,6 +28,7 @@ class TTree;
 class TTreeSRedirector;
 class  AliRawEventHeaderBase;
 class AliTPCCalROC;
+class TVirtualFFT;
 
 class AliTPCclustererMI : public TObject{
 public:
@@ -50,13 +51,13 @@ private:
   Float_t  GetSigmaY2(Int_t iz);
   Float_t  GetSigmaZ2(Int_t iz);
   Float_t  FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz);
-  void AddCluster(AliTPCclusterMI &c);  // add the cluster to the array
+  void AddCluster(AliTPCclusterMI &c, Float_t *matrix, Int_t pos);  // add the cluster to the array
   void UnfoldCluster(Float_t * matrix[7], Float_t recmatrix[5][5], 
                     Float_t & meani, Float_t & meanj, Float_t & sum, Float_t &overlap );
   void FindClusters(AliTPCCalROC * noiseROC);
   Double_t  ProcesSignal(Float_t * signal, Int_t nchannels, Int_t id[3], Double_t &rms, Double_t &pedestalCalib);
-  void DumpHistos();
-
+  void DumpHistos(); 
+  Int_t  TransformFFT(Float_t *input, Float_t threshold, Bool_t locMax, Float_t *freq, Float_t *re, Float_t *im, Float_t *mag, Float_t *phi);
 
   Float_t * fBins;       //!digits array
   Float_t * fResBins;    //!digits array with res. after 1 finder
@@ -86,6 +87,7 @@ private:
   TObjArray * fAmplitudeHisto;          //! array of histograms of amplitudes
   TTreeSRedirector *fDebugStreamer;     //!debug streamer
   const AliTPCRecoParam  * fRecoParam;        //! reconstruction parameters
+  TVirtualFFT *fFFTr2c;                 //! Fast Furier transform object   
   ClassDef(AliTPCclustererMI,1)  // Time Projection Chamber digits
 };