]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCclustererMI.cxx
Adding pad type correction (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.cxx
index 7f4f553de5a42b4c0913e03beb4bc43fab559e73..c932dd4c49bf2a718990b7aa4af3d67fb8951727 100644 (file)
@@ -25,6 +25,7 @@
 //  2. The Output data
 //      2.a TTree with clusters - if  SetOutput(TTree * tree) invoked
 //      2.b TObjArray           - Faster option for HLT
+//      2.c TClonesArray        - Faster option for HLT (smaller memory consumption), activate with fBClonesArray flag
 //
 //  3. Reconstruction setup
 //     see AliTPCRecoParam for list of parameters 
@@ -43,6 +44,7 @@
 #include <TGraph.h>
 #include <TH1F.h>
 #include <TObjArray.h>
+#include <TClonesArray.h>
 #include <TRandom.h>
 #include <TTree.h>
 #include <TTreeStream.h>
@@ -95,13 +97,16 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar
   fInput(0),
   fOutput(0),
   fOutputArray(0),
+  fOutputClonesArray(0),
   fRowCl(0),
   fRowDig(0),
   fParam(0),
   fNcluster(0),
+  fNclusters(0),
   fDebugStreamer(0),
   fRecoParam(0),
-  fBDumpSignal(kFALSE)
+  fBDumpSignal(kFALSE),
+  fBClonesArray(kFALSE)
 {
   //
   // COSNTRUCTOR
@@ -117,8 +122,12 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar
     fRecoParam = AliTPCReconstructor::GetRecoParam();
     if (!fRecoParam)  fRecoParam = AliTPCRecoParam::GetLowFluxParam();
   }
-  fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
-  Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
+  if(AliTPCReconstructor::StreamLevel()>0) {
+    fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
+  }
+
+  //  Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
   fRowCl= new AliTPCClustersRow();
   fRowCl->SetClass("AliTPCclusterMI");
   fRowCl->SetArray(1);
@@ -148,13 +157,16 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
   fInput(0),
   fOutput(0),
   fOutputArray(0),
+  fOutputClonesArray(0),
   fRowCl(0),
   fRowDig(0),
   fParam(0),
   fNcluster(0),
+  fNclusters(0),
   fDebugStreamer(0),
   fRecoParam(0),
-  fBDumpSignal(kFALSE)
+  fBDumpSignal(kFALSE),
+  fBClonesArray(kFALSE)
 {
   //
   // dummy
@@ -180,6 +192,10 @@ AliTPCclustererMI::~AliTPCclustererMI(){
     //fOutputArray->Delete();
     delete fOutputArray;
   }
+  if (fOutputClonesArray){
+    fOutputClonesArray->Delete();
+    delete fOutputClonesArray;
+  }
 }
 
 void AliTPCclustererMI::SetInput(TTree * tree)
@@ -219,10 +235,10 @@ void AliTPCclustererMI::FillRow(){
   //          TObjArray
   //
   if (fOutput) fOutput->Fill();
-  if (!fOutput){
+  if (!fOutput && !fBClonesArray){
     //
     if (!fOutputArray) fOutputArray = new TObjArray(fParam->GetNRowsTotal());
-    if (fRowCl) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID());
+    if (fRowCl && fRowCl->GetArray()->GetEntriesFast()>0) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID());
   }
 }
 
@@ -393,8 +409,6 @@ AliTPCclusterMI &c)
   //unfolding 2
   meani-=i0;
   meanj-=j0;
-  if (gDebug>4)
-    printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
 }
 
 
@@ -522,8 +536,6 @@ void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5]
          }
       }
   }
-  if (gDebug>4) 
-    printf("%f\n", recmatrix[2][2]);
   
 }
 
@@ -591,19 +603,23 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t p
   c.SetZ(x[2]);
   //
   //
-  if (!fRecoParam->GetBYMirror()){
-    if (fSector%36>17){
-      c.SetY(-c.GetY());
-    }
-  }
-
   if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
     c.SetType(-(c.GetType()+3));  //edge clusters
   }
   if (fLoop==2) c.SetType(100);
+  if (!AcceptCluster(&c)) return;
+
+  // select output 
+  TClonesArray * arr = 0;
+  AliTPCclusterMI * cl = 0;
+
+  if(fBClonesArray==kFALSE) {
+     arr = fRowCl->GetArray();
+     cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
+  } else {
+     cl = new ((*fOutputClonesArray)[fNclusters+fNcluster]) AliTPCclusterMI(c);
+  }
 
-  TClonesArray * arr = fRowCl->GetArray();
-  AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
   // if (fRecoParam->DumpSignal() &&matrix ) {
 //     Int_t nbins=0;
 //     Float_t *graph =0;
@@ -617,6 +633,12 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t p
   if (!fRecoParam->DumpSignal()) {
     cl->SetInfo(0);
   }
+  
+  if (AliTPCReconstructor::StreamLevel()>1) {
+     (*fDebugStreamer)<<"Clusters"<<
+       "Cl.="<<cl<<
+       "\n";
+  }
 
   fNcluster++;
 }
@@ -689,7 +711,11 @@ void AliTPCclustererMI::Digits2Clusters()
        Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
        Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
        Int_t bin = i*fMaxTime+j;
-       fBins[bin]=dig/gain;
+       if (gain>0){
+         fBins[bin]=dig/gain;
+       }else{
+         fBins[bin]=0;
+       }
        fSigBins[fNSigBins++]=bin;
       } while (digarr.Next());
     digarr.ExpandTrackBuffer();
@@ -698,6 +724,7 @@ void AliTPCclustererMI::Digits2Clusters()
     FillRow();
     fRowCl->GetArray()->Clear();    
     nclusters+=fNcluster;    
+
     delete[] fBins;
     delete[] fSigBins;
   }  
@@ -729,8 +756,10 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
     fEventType = fEventHeader->Get("Type");  
   }
 
-
-  Int_t nclusters  = 0;
+  // creaate one TClonesArray for all clusters
+  if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
+  // reset counter
+  fNclusters  = 0;
   
   fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
   const Int_t kNIS = fParam->GetNInnerSector();
@@ -762,14 +791,6 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
   //
   for(fSector = 0; fSector < kNS; fSector++) {
 
-    AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
-    AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector);  // pedestal per given sector
-    AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector);  // noise per given sector
-    //check the presence of the calibration
-    if (!noiseROC ||!pedestalROC ) {
-      AliError(Form("Missing calibration per sector\t%d\n",fSector));
-      continue;
-    }
     Int_t nRows = 0;
     Int_t nDDLs = 0, indexDDL = 0;
     if (fSector < kNIS) {
@@ -785,6 +806,23 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
     }
 
+    // load the raw data for corresponding DDLs
+    rawReader->Reset();
+    rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
+
+    // select only good sector 
+    input.Next();
+    if(input.GetSector() != fSector) continue;
+
+    AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
+    AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector);  // pedestal per given sector
+    AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector);  // noise per given sector
+    //check the presence of the calibration
+    if (!noiseROC ||!pedestalROC ) {
+      AliError(Form("Missing calibration per sector\t%d\n",fSector));
+      continue;
+    }
+    
     for (Int_t iRow = 0; iRow < nRows; iRow++) {
       Int_t maxPad;
       if (fSector < kNIS)
@@ -797,20 +835,23 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
       allNSigBins[iRow] = 0;
     }
     
-    // Loas the raw data for corresponding DDLs
-    rawReader->Reset();
-    rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
     Int_t digCounter=0;
     // Begin loop over altro data
     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
     Float_t gain =1;
     Int_t lastPad=-1;
+
+    input.Reset();
     while (input.Next()) {
       if (input.GetSector() != fSector)
        AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
 
       
       Int_t iRow = input.GetRow();
+      if (iRow < 0){
+       continue;
+      }
+
       if (iRow < 0 || iRow >= nRows){
        AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
                      iRow, 0, nRows -1));
@@ -840,9 +881,14 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
       //signal
       Float_t signal = input.GetSignal();
       if (!calcPedestal && signal <= zeroSup) continue;      
+
       if (!calcPedestal) {
        Int_t bin = iPad*fMaxTime+iTimeBin;
-       allBins[iRow][bin] = signal/gain;
+       if (gain>0){
+         allBins[iRow][bin] = signal/gain;
+       }else{
+         allBins[iRow][bin] =0;
+       }
        allSigBins[iRow][allNSigBins[iRow]++] = bin;
       }else{
        allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
@@ -907,6 +953,48 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
        }
       }
     }
+
+    if (AliTPCReconstructor::StreamLevel()>3) {
+      for (Int_t iRow = 0; iRow < nRows; iRow++) {
+       Int_t maxPad;
+       if (fSector < kNIS)
+         maxPad = fParam->GetNPadsLow(iRow);
+       else
+         maxPad = fParam->GetNPadsUp(iRow);
+       
+       for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
+         for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
+           Int_t bin = iPad*fMaxTime+iTimeBin;
+           Float_t signal = allBins[iRow][bin];
+           if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
+             Double_t x[]={iRow,iPad-3,iTimeBin-3};
+             Int_t i[]={fSector};
+             AliTPCTransform trafo;
+             trafo.Transform(x,i,0,1);
+             Double_t gx[3]={x[0],x[1],x[2]};
+             trafo.RotatedGlobal2Global(fSector,gx);
+             
+              if (AliTPCReconstructor::StreamLevel()>0) {
+             (*fDebugStreamer)<<"Digits"<<
+               "sec="<<fSector<<
+               "row="<<iRow<<
+               "pad="<<iPad<<
+               "time="<<iTimeBin<<
+               "sig="<<signal<<
+               "x="<<x[0]<<
+               "y="<<x[1]<<
+               "z="<<x[2]<<      
+               "gx="<<gx[0]<<
+               "gy="<<gx[1]<<
+               "gz="<<gx[2]<<
+               "\n";
+             }
+           }
+         }
+       }
+      }
+    }
+
     // Now loop over rows and find clusters
     for (fRow = 0; fRow < nRows; fRow++) {
       fRowCl->SetID(fParam->GetIndex(fSector, fRow));
@@ -926,12 +1014,14 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
       fNSigBins = allNSigBins[fRow];
 
       FindClusters(noiseROC);
+
       FillRow();
-      fRowCl->GetArray()->Clear();    
-      nclusters += fNcluster;    
+      if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear();    
+      fNclusters += fNcluster;    
+
     } // End of loop to find clusters
   } // End of loop over sectors
-  
+
   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
     delete [] allBins[iRow];
     delete [] allSigBins[iRow];
@@ -940,13 +1030,17 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
   delete [] allSigBins;
   delete [] allNSigBins;
   
-//   if (rawReader->GetEventId() && fOutput ){
-//     Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
-//   }else{
-//     Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), nclusters);
-    
-//   }
+  if (rawReader->GetEventId() && fOutput ){
+    Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
+  }
   
+  if(rawReader->GetEventId()) {
+    Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
+  }
+
+  if(fBClonesArray) {
+    //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
+  }
 }
 
 void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
@@ -1001,6 +1095,22 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
   }
 }
 
+Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
+  //
+  // Currently hack to filter digital noise (15.06.2008)
+  // To be parameterized in the AliTPCrecoParam
+  // More inteligent way  to be used in future
+  // Acces to the proper pedestal file needed
+  //
+  if (cl->GetMax()<400) return kTRUE;
+  Double_t ratio = cl->GetQ()/cl->GetMax();
+  if (cl->GetMax()>700){
+    if ((ratio - int(ratio)>0.8)) return kFALSE;
+  }
+  if ((ratio - int(ratio)<0.95)) return kTRUE;
+  return kFALSE;
+}
+
 
 Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
   //
@@ -1014,7 +1124,6 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
   // 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;
@@ -1025,7 +1134,8 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
   Int_t    firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
   //
   UShort_t histo[kPedMax];
-  memset(histo,0,kPedMax*sizeof(UShort_t));
+  //memset(histo,0,kPedMax*sizeof(UShort_t));
+  for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
   for (Int_t i=0; i<fMaxTime; i++){
     if (signal[i]<=0) continue;
     if (signal[i]>max && i>firstBin) {
@@ -1096,7 +1206,8 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
   //
   // Dump mean signal info
   //
-  (*fDebugStreamer)<<"Signal"<<
+    if (AliTPCReconstructor::StreamLevel()>0) {
+    (*fDebugStreamer)<<"Signal"<<
     "TimeStamp="<<fTimeStamp<<
     "EventType="<<fEventType<<
     "Sector="<<uid[0]<<
@@ -1115,11 +1226,10 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
     "RMSCalib="<<rmsCalib<<
     "PedCalib="<<pedestalCalib<<
     "\n";
+    }
   //
   // fill pedestal histogram
   //
-  AliTPCROC * roc = AliTPCROC::Instance();
-
   //
   //
   //
@@ -1131,10 +1241,11 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
     dsignal[i] = signal[i];
   }
 
-  TGraph * graph;
+  TGraph * graph=0;
   //
   // Big signals dumping
   //    
+  if (AliTPCReconstructor::StreamLevel()>0) {
   if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin()) 
     (*fDebugStreamer)<<"SignalB"<<     // pads with signal
       "TimeStamp="<<fTimeStamp<<
@@ -1155,6 +1266,7 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
       "RMS09="<<rms09<<
       "\n";
   delete graph;  
+  }
 
   delete [] dsignal;
   delete [] dtime;