From 3e1a3ad8f2e5e0e979c9eec471ff0a374e5073d5 Mon Sep 17 00:00:00 2001 From: cblume Date: Mon, 7 May 2001 08:06:44 +0000 Subject: [PATCH] Speedup of the code. Create only AliTRDcluster --- TRD/AliTRDclusterizer.cxx | 193 ++++++++++++----- TRD/AliTRDclusterizer.h | 14 +- TRD/AliTRDclusterizerV0.cxx | 17 +- TRD/AliTRDclusterizerV1.cxx | 403 ++++++++++++++++-------------------- TRD/AliTRDclusterizerV1.h | 17 +- 5 files changed, 339 insertions(+), 305 deletions(-) diff --git a/TRD/AliTRDclusterizer.cxx b/TRD/AliTRDclusterizer.cxx index 8c6b7ef8e8d..f09587ec489 100644 --- a/TRD/AliTRDclusterizer.cxx +++ b/TRD/AliTRDclusterizer.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.7 2001/03/30 14:40:14 cblume +Update of the digitization parameter + Revision 1.6 2000/11/01 14:53:20 cblume Merge with TRD-develop @@ -76,6 +79,7 @@ Add new TRD classes #include "AliRun.h" #include "AliTRD.h" #include "AliTRDclusterizer.h" +#include "AliTRDcluster.h" #include "AliTRDrecPoint.h" #include "AliTRDgeometry.h" @@ -88,8 +92,11 @@ AliTRDclusterizer::AliTRDclusterizer():TNamed() // AliTRDclusterizer default constructor // - fInputFile = NULL; - fEvent = 0; + fInputFile = NULL; + fOutputFile = NULL; + fClusterTree = NULL; + fEvent = 0; + fVerbose = 0; } @@ -101,10 +108,11 @@ AliTRDclusterizer::AliTRDclusterizer(const Text_t* name, const Text_t* title) // AliTRDclusterizer default constructor // - fInputFile = NULL; - fEvent = 0; - - Init(); + fInputFile = NULL; + fOutputFile = NULL; + fClusterTree = NULL; + fEvent = 0; + fVerbose = 0; } @@ -131,6 +139,15 @@ AliTRDclusterizer::~AliTRDclusterizer() delete fInputFile; } + if (fOutputFile) { + fOutputFile->Close(); + delete fOutputFile; + } + + if (fClusterTree) { + delete fClusterTree; + } + } //_____________________________________________________________________________ @@ -152,22 +169,76 @@ void AliTRDclusterizer::Copy(TObject &c) // Copy function // - ((AliTRDclusterizer &) c).fInputFile = NULL; - ((AliTRDclusterizer &) c).fEvent = 0; + ((AliTRDclusterizer &) c).fInputFile = NULL; + ((AliTRDclusterizer &) c).fOutputFile = NULL; + ((AliTRDclusterizer &) c).fClusterTree = NULL; + ((AliTRDclusterizer &) c).fEvent = 0; + ((AliTRDclusterizer &) c).fVerbose = fVerbose; } //_____________________________________________________________________________ -void AliTRDclusterizer::Init() +Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent) { // - // Initializes the cluster finder + // Opens the AliROOT file. Output and input are in the same file // + OpenInput(name,nEvent); + OpenOutput(name); + + return kTRUE; + } //_____________________________________________________________________________ -Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent) +Bool_t AliTRDclusterizer::Open(const Char_t *inname, const Char_t *outname + , Int_t nEvent) +{ + // + // Opens the AliROOT file. Output and input are in different files + // + + OpenInput(inname,nEvent); + OpenOutput(outname); + + return kTRUE; + +} + +//_____________________________________________________________________________ +Bool_t AliTRDclusterizer::OpenOutput(const Char_t *name) +{ + // + // Open the output file + // + + TDirectory *savedir = NULL; + + if (!fInputFile) return kFALSE; + + if (strcmp(name,fInputFile->GetName()) != 0) { + savedir = gDirectory; + printf("AliTRDclusterizer::OpenOutput -- "); + printf("Open the output file %s.\n",name); + fOutputFile = new TFile(name,"RECREATE"); + } + + // Create a tree for the cluster + fClusterTree = new TTree("ClusterTree","Tree with TRDcluster"); + TObjArray *ioArray = 0; + fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0); + + if (savedir) { + savedir->cd(); + } + + return kTRUE; + +} + +//_____________________________________________________________________________ +Bool_t AliTRDclusterizer::OpenInput(const Char_t *name, Int_t nEvent) { // // Opens a ROOT-file with TRD-hits and reads in the digits-tree @@ -176,23 +247,21 @@ Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent) // Connect the AliRoot file containing Geometry, Kine, and Hits fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(name); if (!fInputFile) { - printf("AliTRDclusterizer::Open -- "); + printf("AliTRDclusterizer::OpenInput -- "); printf("Open the ALIROOT-file %s.\n",name); fInputFile = new TFile(name,"UPDATE"); } else { - printf("AliTRDclusterizer::Open -- "); + printf("AliTRDclusterizer::OpenInput -- "); printf("%s is already open.\n",name); } - // Get AliRun object from file or create it if not on file + // Get AliRun object from file + gAlice = (AliRun *) fInputFile->Get("gAlice"); if (!(gAlice)) { - gAlice = (AliRun *) fInputFile->Get("gAlice"); - if (!(gAlice)) { - printf("AliTRDclusterizer::Open -- "); - printf("Could not find AliRun object.\n"); - return kFALSE; - } + printf("AliTRDclusterizer::OpenInput -- "); + printf("Could not find AliRun object.\n"); + return kFALSE; } fEvent = nEvent; @@ -200,17 +269,18 @@ Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent) // Import the Trees for the event nEvent in the file Int_t nparticles = gAlice->GetEvent(fEvent); if (nparticles <= 0) { - printf("AliTRDclusterizer::Open -- "); + printf("AliTRDclusterizer::OpenInput -- "); printf("No entries in the trees for event %d.\n",fEvent); return kFALSE; } - // Create a tree for the reconstructed points - TTree *recPointTree = new TTree("ClusterTree","Tree with clusters and rec. points"); - TObjArray *ioArray = 0; - recPointTree->Branch("TRDrecPoints","TObjArray",&ioArray,32000,0); -// TObjArray *iopointer = 0; -// recPointTree->Branch("Clusters","TObjArray",&iopointer,32000,0); + // Get the TRD object + fTRD = (AliTRD*) gAlice->GetDetector("TRD"); + if (!fTRD) { + printf("AliTRDclusterizer::OpenInput -- "); + printf("No TRD detector object found\n"); + return kFALSE; + } return kTRUE; @@ -220,59 +290,68 @@ Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent) Bool_t AliTRDclusterizer::WriteClusters(Int_t det) { // - // Fills TRDrecPoints branch in TRDrecPoints## tree with rec. points + // Fills TRDcluster branch in the tree with the clusters // found in detector = det. For det=-1 writes the tree. - // For det=-2 recreates the tree. // - Char_t treeName[14]; - sprintf(treeName,"TRDrecPoints%d", fEvent); - - if (det == -2) { - fInputFile->Delete(treeName); - TTree *tree = new TTree(treeName,"Tree with TRD rec. points"); - tree->Write(); - return kTRUE; + if ((det < -1) || (det >= AliTRDgeometry::Ndet())) { + printf("AliTRDclusterizer::WriteClusters -- "); + printf("Unexpected detector index %d.\n",det); + return kFALSE; } + + TDirectory *savedir = gDirectory; - //TTree *tree = (TTree *) fInputFile->Get(treeName); - TTree *tree = (TTree *) fInputFile->Get("ClusterTree"); - TBranch *branch = tree->GetBranch("TRDrecPoints"); + if (fOutputFile) { + fOutputFile->cd(); + } - if(!branch) { + TBranch *branch = fClusterTree->GetBranch("TRDcluster"); + if (!branch) { TObjArray *ioArray = 0; - branch = tree->Branch("TRDrecPoints","TObjArray",&ioArray,32000,0); + branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0); } if ((det >= 0) && (det < AliTRDgeometry::Ndet())) { - AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD"); - Int_t nRecPoints = TRD->RecPoints()->GetEntriesFast(); - TObjArray *fDetRecPoints = new TObjArray(400); - - for (Int_t i=0; iRecPoints()->UncheckedAt(i); - if(det == p->GetDetector()) fDetRecPoints->AddLast(p); - else printf("attempt to write a RecPoint with unexpected detector index"); + Int_t nRecPoints = fTRD->RecPoints()->GetEntriesFast(); + TObjArray *detRecPoints = new TObjArray(400); + + for (Int_t i = 0; i < nRecPoints; i++) { + AliTRDcluster *c = (AliTRDcluster *) fTRD->RecPoints()->UncheckedAt(i); + if (det == c->GetDetector()) { + detRecPoints->AddLast(c); + } + else { + printf("AliTRDclusterizer::WriteClusters --"); + printf("Attempt to write a cluster with unexpected detector index\n"); + } } - branch->SetAddress(&fDetRecPoints); - tree->Fill(); + branch->SetAddress(&detRecPoints); + fClusterTree->Fill(); + return kTRUE; + } if (det == -1) { - printf("\rAliTRDclusterizer::WriteClusters -- "); + printf("AliTRDclusterizer::WriteClusters -- "); printf("Writing the cluster tree %-18s for event %d.\n" - ,tree->GetName(),fEvent); + ,fClusterTree->GetName(),fEvent); - tree->Write(); + fClusterTree->Write(); + return kTRUE; + } - printf("\rAliTRDclusterizer::WriteClusters -- "); - printf("Unexpected detector index %d.\n", det); + savedir->cd(); + + printf("AliTRDclusterizer::WriteClusters -- "); + printf("Unexpected detector index %d.\n",det); + return kFALSE; } diff --git a/TRD/AliTRDclusterizer.h b/TRD/AliTRDclusterizer.h index af8997d4583..94a1532d19a 100644 --- a/TRD/AliTRDclusterizer.h +++ b/TRD/AliTRDclusterizer.h @@ -24,18 +24,26 @@ class AliTRDclusterizer : public TNamed { AliTRDclusterizer &operator=(const AliTRDclusterizer &c); virtual void Copy(TObject &c); - virtual void Init(); virtual Bool_t Open(const Char_t *name, Int_t nEvent = 0); + virtual Bool_t Open(const Char_t *inname, const Char_t *outname, Int_t nEvent = 0); + virtual Bool_t OpenInput(const Char_t *name, Int_t nEvent = 0); + virtual Bool_t OpenOutput(const Char_t *name); virtual Bool_t MakeClusters() = 0; virtual Bool_t WriteClusters(Int_t det); + void SetVerbose(Int_t v = 1) { fVerbose = v; }; + protected: TFile *fInputFile; //! AliROOT input file - + TFile *fOutputFile; //! AliROOT output file + TTree *fClusterTree; //! Tree with the cluster + AliTRD *fTRD; //! The TRD object + Int_t fEvent; //! Event number + Int_t fVerbose; // Sets the verbose level - ClassDef(AliTRDclusterizer,1) // TRD-Cluster manager base class + ClassDef(AliTRDclusterizer,2) // TRD-Cluster manager base class }; diff --git a/TRD/AliTRDclusterizerV0.cxx b/TRD/AliTRDclusterizerV0.cxx index 6162866d134..535bf4ae2aa 100644 --- a/TRD/AliTRDclusterizerV0.cxx +++ b/TRD/AliTRDclusterizerV0.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.7 2001/02/14 18:22:26 cblume +Change in the geometry of the padplane + Revision 1.6 2000/11/01 14:53:20 cblume Merge with TRD-develop @@ -143,16 +146,14 @@ Bool_t AliTRDclusterizerV0::MakeClusters() // Generates the cluster // - // Get the pointer to the detector class and check for version 1 - AliTRD *trd = (AliTRD*) gAlice->GetDetector("TRD"); - if (trd->IsVersion() != 0) { + if (fTRD->IsVersion() != 0) { printf("AliTRDclusterizerV0::MakeCluster -- "); printf("TRD must be version 0 (fast simulator).\n"); return kFALSE; } // Get the geometry - AliTRDgeometry *geo = trd->GetGeometry(); + AliTRDgeometry *geo = fTRD->GetGeometry(); printf("AliTRDclusterizerV0::MakeCluster -- "); printf("Start creating cluster.\n"); @@ -193,12 +194,12 @@ Bool_t AliTRDclusterizerV0::MakeClusters() nBytes += hitTree->GetEvent(iTrack); // Get the number of hits in the TRD created by this particle - Int_t nHit = trd->Hits()->GetEntriesFast(); + Int_t nHit = fTRD->Hits()->GetEntriesFast(); // Loop through the TRD hits for (Int_t iHit = 0; iHit < nHit; iHit++) { - if (!(hit = (AliTRDhit *) trd->Hits()->UncheckedAt(iHit))) + if (!(hit = (AliTRDhit *) fTRD->Hits()->UncheckedAt(iHit))) continue; Float_t pos[3]; @@ -316,7 +317,7 @@ Bool_t AliTRDclusterizerV0::MakeClusters() Int_t detector = recPoint1->GetDetector(); Int_t digits[3] = {0}; Int_t tr[9] = {-1}; - trd->AddRecPoint(smear,digits,detector,0.0,tr); + fTRD->AddCluster(smear,digits,detector,0.0,tr,0); } @@ -328,7 +329,7 @@ Bool_t AliTRDclusterizerV0::MakeClusters() } printf("AliTRDclusterizerV0::MakeCluster -- "); - printf("Found %d points.\n",trd->RecPoints()->GetEntries()); + printf("Found %d points.\n",fTRD->RecPoints()->GetEntries()); printf("AliTRDclusterizerV0::MakeCluster -- "); printf("Fill the cluster tree.\n"); clusterTree->Fill(); diff --git a/TRD/AliTRDclusterizerV1.cxx b/TRD/AliTRDclusterizerV1.cxx index e76f11e6b70..7a7a279939d 100644 --- a/TRD/AliTRDclusterizerV1.cxx +++ b/TRD/AliTRDclusterizerV1.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.9 2000/11/01 14:53:20 cblume +Merge with TRD-develop + Revision 1.1.4.5 2000/10/15 23:40:01 cblume Remove AliTRDconst @@ -70,7 +73,6 @@ Add new TRD classes #include "AliTRDmatrix.h" #include "AliTRDgeometry.h" #include "AliTRDdigitizer.h" -#include "AliTRDrecPoint.h" #include "AliTRDdataArrayF.h" #include "AliTRDdataArrayI.h" #include "AliTRDdigitsManager.h" @@ -86,6 +88,9 @@ AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer() fDigitsManager = NULL; + fClusMaxThresh = 0; + fClusSigThresh = 0; + } //_____________________________________________________________________________ @@ -147,7 +152,6 @@ void AliTRDclusterizerV1::Copy(TObject &c) ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh; ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh; - ((AliTRDclusterizerV1 &) c).fClusMethod = fClusMethod; ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL; AliTRDclusterizer::Copy(c); @@ -162,9 +166,8 @@ void AliTRDclusterizerV1::Init() // // The default parameter for the clustering - fClusMaxThresh = 5.0; - fClusSigThresh = 2.0; - fClusMethod = 1; + fClusMaxThresh = 3; + fClusSigThresh = 1; } @@ -195,16 +198,14 @@ Bool_t AliTRDclusterizerV1::MakeClusters() Int_t row, col, time; - // Get the pointer to the detector class and check for version 1 - AliTRD *trd = (AliTRD*) gAlice->GetDetector("TRD"); - if (trd->IsVersion() != 1) { + if (fTRD->IsVersion() != 1) { printf("AliTRDclusterizerV1::MakeCluster -- "); printf("TRD must be version 1 (slow simulator).\n"); return kFALSE; } // Get the geometry - AliTRDgeometry *geo = trd->GetGeometry(); + AliTRDgeometry *geo = fTRD->GetGeometry(); printf("AliTRDclusterizerV1::MakeCluster -- "); printf("Start creating clusters.\n"); @@ -214,40 +215,47 @@ Bool_t AliTRDclusterizerV1::MakeClusters() AliTRDdataArrayI *track1; AliTRDdataArrayI *track2; - // Parameters - Float_t maxThresh = fClusMaxThresh; // threshold value for maximum - Float_t signalThresh = fClusSigThresh; // threshold value for digit signal - Int_t clusteringMethod = fClusMethod; // clustering method option (for testing) + // Threshold value for the maximum + Int_t maxThresh = fClusMaxThresh; + // Threshold value for the digit signal + Int_t sigThresh = fClusSigThresh; // Iteration limit for unfolding procedure const Float_t kEpsilon = 0.01; const Int_t kNclus = 3; const Int_t kNsig = 5; + const Int_t kNtrack = 3 * kNclus; + + Float_t padSignal[kNsig]; + Float_t clusterSignal[kNclus]; + Float_t clusterPads[kNclus]; + Int_t clusterDigit[kNclus]; + Int_t clusterTracks[kNtrack]; Int_t chamBeg = 0; Int_t chamEnd = AliTRDgeometry::Ncham(); - if (trd->GetSensChamber() >= 0) { - chamBeg = trd->GetSensChamber(); + if (fTRD->GetSensChamber() >= 0) { + chamBeg = fTRD->GetSensChamber(); chamEnd = chamBeg + 1; } Int_t planBeg = 0; Int_t planEnd = AliTRDgeometry::Nplan(); - if (trd->GetSensPlane() >= 0) { - planBeg = trd->GetSensPlane(); + if (fTRD->GetSensPlane() >= 0) { + planBeg = fTRD->GetSensPlane(); planEnd = planBeg + 1; } Int_t sectBeg = 0; Int_t sectEnd = AliTRDgeometry::Nsect(); - // *** Start clustering *** in every chamber + // Start clustering in every chamber for (Int_t icham = chamBeg; icham < chamEnd; icham++) { for (Int_t iplan = planBeg; iplan < planEnd; iplan++) { for (Int_t isect = sectBeg; isect < sectEnd; isect++) { - if (trd->GetSensSector() >= 0) { - Int_t sens1 = trd->GetSensSector(); - Int_t sens2 = sens1 + trd->GetSensSectorRange(); + if (fTRD->GetSensSector() >= 0) { + Int_t sens1 = fTRD->GetSensSector(); + Int_t sens2 = sens1 + fTRD->GetSensSectorRange(); sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) * AliTRDgeometry::Nsect(); if (sens1 < sens2) { @@ -260,215 +268,155 @@ Bool_t AliTRDclusterizerV1::MakeClusters() Int_t idet = geo->GetDetector(iplan,icham,isect); - Int_t nClusters = 0; + Int_t nClusters = 0; + Int_t nClustersUnfold = 0; + printf("AliTRDclusterizerV1::MakeCluster -- "); printf("Analyzing chamber %d, plane %d, sector %d.\n" - ,icham,iplan,isect); + ,icham,iplan,isect); - Int_t nRowMax = geo->GetRowMax(iplan,icham,isect); - Int_t nColMax = geo->GetColMax(iplan); - Int_t nTimeMax = geo->GetTimeMax(); + Int_t nRowMax = geo->GetRowMax(iplan,icham,isect); + Int_t nColMax = geo->GetColMax(iplan); + Int_t nTimeBefore = geo->GetTimeBefore(); + Int_t nTimeTotal = geo->GetTimeTotal(); - // Create a detector matrix to keep maxima - AliTRDmatrix *digitMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax - ,isect,icham,iplan); - // Create a matrix to contain maximum flags - AliTRDmatrix *maximaMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax - ,isect,icham,iplan); - - // Create a matrix for track indexes - AliTRDmatrix *trackMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax - ,isect,icham,iplan); - - // Read in the digits + // Get the digits digits = fDigitsManager->GetDigits(idet); + digits->Expand(); track0 = fDigitsManager->GetDictionary(idet,0); + track0->Expand(); track1 = fDigitsManager->GetDictionary(idet,1); + track1->Expand(); track2 = fDigitsManager->GetDictionary(idet,2); - - // Loop through the detector pixel - for (time = 0; time < nTimeMax; time++) { - for ( col = 0; col < nColMax; col++) { - for ( row = 0; row < nRowMax; row++) { - - Int_t signal = digits->GetData(row,col,time); - Int_t index = digits->GetIndex(row,col,time); - Int_t t[3] = {-1}; - t[0] = track0->GetData(row,col,time) - 1; - t[1] = track1->GetData(row,col,time) - 1; - t[2] = track2->GetData(row,col,time) - 1; - - // Fill the detector matrix - if (signal > signalThresh) { - // Store the signal amplitude - digitMatrix->SetSignal(row,col,time,signal); - // Store the digits number - digitMatrix->AddTrack(row,col,time,index); - for(Int_t i = 0; i < 3; i++) { - trackMatrix->AddTrack(row,col,time,t[i]); - } - } - } - } - } - - // Loop chamber and find maxima in digitMatrix - for ( row = 0; row < nRowMax; row++) { - for ( col = 1; col < nColMax; col++) { - for (time = 0; time < nTimeMax; time++) { - - if (digitMatrix->GetSignal(row,col,time) - < digitMatrix->GetSignal(row,col - 1,time)) { - // really maximum? - if (col > 1) { - if (digitMatrix->GetSignal(row,col - 2,time) - < digitMatrix->GetSignal(row,col - 1,time)) { - // yes, so set maximum flag - maximaMatrix->SetSignal(row,col - 1,time,1); - } - else maximaMatrix->SetSignal(row,col - 1,time,0); - } - } - - } // time - } // col - } // row - - // now check maxima and calculate cluster position - for ( row = 0; row < nRowMax; row++) { - for ( col = 1; col < nColMax; col++) { - for (time = 0; time < nTimeMax; time++) { - - if ((maximaMatrix->GetSignal(row,col,time) > 0) - && (digitMatrix->GetSignal(row,col,time) > maxThresh)) { - - // Ratio resulting from unfolding - Float_t ratio = 0; - // Signals on max and neighbouring pads - Float_t padSignal[kNsig] = {0}; - // Signals from cluster - Float_t clusterSignal[kNclus] = {0}; - // Cluster pad info - Float_t clusterPads[kNclus] = {0}; - // Cluster digit info - Int_t clusterDigit[kNclus] = {0}; - // Cluster MC tracks info - const Int_t nt = kNclus*3; - Int_t clusterTracks[nt] = {-1}; + track2->Expand(); + + // Loop through the chamber and find the maxima + for ( row = 0; row < nRowMax; row++) { + for ( col = 2; col < nColMax; col++) { + for (time = 0; time < nTimeTotal; time++) { + + Int_t signalL = digits->GetDataUnchecked(row,col ,time); + Int_t signalM = digits->GetDataUnchecked(row,col-1,time); + Int_t signalR = digits->GetDataUnchecked(row,col-2,time); + + // Look for the maximum + if ((signalM >= maxThresh) && + (signalL >= sigThresh) && + (signalR >= sigThresh)) { + if ((signalL < signalM) && (signalR < signalM)) { + // Maximum found, mark the position by a negative signal + digits->SetDataUnchecked(row,col-1,time,-signalM); + } + } + + } + } + } + + // Now check the maxima and calculate the cluster position + for ( row = 0; row < nRowMax ; row++) { + for ( col = 1; col < nColMax-1; col++) { + for (time = 0; time < nTimeTotal; time++) { + + // Maximum found ? + if (digits->GetDataUnchecked(row,col,time) < 0) { Int_t iPad; for (iPad = 0; iPad < kNclus; iPad++) { - clusterSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time); - clusterDigit[iPad] = digitMatrix->GetTrack(row,col-1+iPad,time,0); - for (Int_t j = 0; j < 3; j++) { - clusterTracks[iPad*3+j] = trackMatrix->GetTrack(row,col-1+iPad,time,j); - } + Int_t iPadCol = col - 1 + iPad; + clusterSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row + ,iPadCol + ,time)); + clusterDigit[iPad] = digits->GetIndexUnchecked(row,iPadCol,time); + clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1; + clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1; + clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1; } - // neighbouring maximum on right side? - if (col < nColMax - 2) { - if (maximaMatrix->GetSignal(row,col + 2,time) > 0) { - - for (iPad = 0; iPad < 5; iPad++) { - padSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time); + // Neighbouring maximum on right side? + Int_t iType = 0; + if (col < nColMax-3) { + if (digits->GetDataUnchecked(row,col+2,time) < 0) { + for (iPad = 0; iPad < kNsig; iPad++) { + padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row + ,col-1+iPad + ,time)); } - - // unfold: - ratio = Unfold(kEpsilon, padSignal); - - // set signal on overlapping pad to ratio - clusterSignal[2] *= ratio; - + // Unfold the two maxima and set the signal on + // the overlapping pad to the ratio + clusterSignal[2] *= Unfold(kEpsilon,padSignal); + iType = 1; } } - - // Calculate the position of the cluster - switch (clusteringMethod) { - case 1: - // method 1: simply center of mass - clusterPads[0] = row + 0.5; - clusterPads[1] = col + 0.5 + (clusterSignal[2] - clusterSignal[0]) / - (clusterSignal[0] + clusterSignal[1] + clusterSignal[2]); - clusterPads[2] = time + 0.5; - - -/* printf("col = %d, left = %f, center = %f, right = %f, - final =%f \n", col, - digitMatrix->GetSignal(row,col-1,time), - digitMatrix->GetSignal(row,col,time), - digitMatrix->GetSignal(row,col+1,time), - clusterPads[1]); - - printf("col = %d, sig(0) = %f, sig(1) = %f, sig(2) = %f, - final =%f \n", col, - clusterSignal[0], - clusterSignal[1], - clusterSignal[2], - clusterPads[1]); - -*/ - - nClusters++; - break; - case 2: - // method 2: integral gauss fit on 3 pads - TH1F *hPadCharges = new TH1F("hPadCharges", "Charges on center 3 pads" - , 5, -1.5, 3.5); - for (Int_t iCol = -1; iCol <= 3; iCol++) { - if (clusterSignal[iCol] < 1) clusterSignal[iCol] = 1; - hPadCharges->Fill(iCol, clusterSignal[iCol]); - } - hPadCharges->Fit("gaus", "IQ", "SAME", -0.5, 2.5); - TF1 *fPadChargeFit = hPadCharges->GetFunction("gaus"); - Double_t colMean = fPadChargeFit->GetParameter(1); - clusterPads[0] = row + 0.5; - clusterPads[1] = col - 1.5 + colMean; - clusterPads[2] = time + 0.5; - - delete hPadCharges; - - nClusters++; - break; + Float_t clusterCharge = clusterSignal[0] + + clusterSignal[1] + + clusterSignal[2]; + + // Calculate the position of the cluster by using the + // center of gravity method + clusterPads[0] = row + 0.5; + clusterPads[1] = col + 0.5 + + (clusterSignal[2] - clusterSignal[0]) + / clusterCharge; + // Take the shift of the additional time bins into account + clusterPads[2] = time - nTimeBefore + 0.5; + + if (fVerbose) { + printf("-----------------------------------------------------------\n"); + printf("Create cluster no. %d\n",nClusters); + printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0] + ,clusterPads[1] + ,clusterPads[2]); + printf("Indices: %d, %d, %d\n",clusterDigit[0] + ,clusterDigit[1] + ,clusterDigit[2]); + printf("Total charge = %f\n",clusterCharge); + printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0] + ,clusterTracks[1] + ,clusterTracks[2]); + printf(" pad1 %d, %d, %d\n",clusterTracks[3] + ,clusterTracks[4] + ,clusterTracks[5]); + printf(" pad2 %d, %d, %d\n",clusterTracks[6] + ,clusterTracks[7] + ,clusterTracks[8]); } - Float_t clusterCharge = clusterSignal[0] - + clusterSignal[1] - + clusterSignal[2]; + nClusters++; + if (iType == 1) nClustersUnfold++; // Add the cluster to the output array - trd->AddRecPoint(clusterPads,clusterDigit,idet,clusterCharge,clusterTracks); + fTRD->AddCluster(clusterPads + ,clusterDigit + ,idet + ,clusterCharge + ,clusterTracks + ,iType); } - } // time - } // col - } // row + } + } + } - printf("AliTRDclusterizerV1::MakeCluster -- "); - printf("Number of clusters found: %d\n",nClusters); + // Compress the arrays + digits->Compress(1,0); + track0->Compress(1,0); + track1->Compress(1,0); + track2->Compress(1,0); + // Write the cluster and reset the array WriteClusters(idet); - trd->ResetRecPoints(); + fTRD->ResetRecPoints(); - delete digitMatrix; - delete maximaMatrix; - delete trackMatrix; - - } // isect - } // iplan - } // icham - -// printf("AliTRDclusterizerV1::MakeCluster -- "); -// printf("Total number of points found: %d\n" -// ,trd->RecPoints()->GetEntries()); + printf("AliTRDclusterizerV1::MakeCluster -- "); + printf("Found %d (%d unfolded) clusters.\n" + ,nClusters,nClustersUnfold); -// // Get the pointer to the cluster branch -// TTree *clusterTree = gAlice->TreeR(); + } + } + } -// // Fill the cluster-branch -// printf("AliTRDclusterizerV1::MakeCluster -- "); -// printf("Fill the cluster tree.\n"); -// clusterTree->Fill(); printf("AliTRDclusterizerV1::MakeCluster -- "); printf("Done.\n"); @@ -486,41 +434,41 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal) // The resulting ratio is then returned to the calling method. // - Int_t itStep = 0; // count iteration steps + Int_t itStep = 0; // Count iteration steps - Float_t ratio = 0.5; // start value for ratio - Float_t prevRatio = 0; // store previous ratio + Float_t ratio = 0.5; // Start value for ratio + Float_t prevRatio = 0; // Store previous ratio - Float_t newLeftSignal[3] = {0}; // array to store left cluster signal - Float_t newRightSignal[3] = {0}; // array to store right cluster signal + Float_t newLeftSignal[3] = {0}; // Array to store left cluster signal + Float_t newRightSignal[3] = {0}; // Array to store right cluster signal - // start iteration: + // Start the iteration while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) { itStep++; prevRatio = ratio; - // cluster position according to charge ratio - Float_t maxLeft = (ratio*padSignal[2] - padSignal[0]) / - (padSignal[0] + padSignal[1] + ratio*padSignal[2]); - Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) / - ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]); + // Cluster position according to charge ratio + Float_t maxLeft = (ratio*padSignal[2] - padSignal[0]) + / (padSignal[0] + padSignal[1] + ratio*padSignal[2]); + Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) + / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]); - // set cluster charge ratio + // Set cluster charge ratio Float_t ampLeft = padSignal[1]; Float_t ampRight = padSignal[3]; - // apply pad response to parameters - newLeftSignal[0] = ampLeft*PadResponse(-1 - maxLeft); - newLeftSignal[1] = ampLeft*PadResponse( 0 - maxLeft); - newLeftSignal[2] = ampLeft*PadResponse( 1 - maxLeft); + // Apply pad response to parameters + newLeftSignal[0] = ampLeft * PadResponse(-1 - maxLeft); + newLeftSignal[1] = ampLeft * PadResponse( 0 - maxLeft); + newLeftSignal[2] = ampLeft * PadResponse( 1 - maxLeft); - newRightSignal[0] = ampRight*PadResponse(-1 - maxRight); - newRightSignal[1] = ampRight*PadResponse( 0 - maxRight); - newRightSignal[2] = ampRight*PadResponse( 1 - maxRight); + newRightSignal[0] = ampRight * PadResponse(-1 - maxRight); + newRightSignal[1] = ampRight * PadResponse( 0 - maxRight); + newRightSignal[2] = ampRight * PadResponse( 1 - maxRight); - // calculate new overlapping ratio - ratio = newLeftSignal[2]/(newLeftSignal[2] + newRightSignal[0]); + // Calculate new overlapping ratio + ratio = newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]); } @@ -535,15 +483,16 @@ Float_t AliTRDclusterizerV1::PadResponse(Float_t x) // The pad response for the chevron pads. // We use a simple Gaussian approximation which should be good // enough for our purpose. + // Updated for new PRF 1/5/01. // // The parameters for the response function - const Float_t kA = 0.8872; - const Float_t kB = -0.00573; - const Float_t kC = 0.454; - const Float_t kC2 = kC*kC; + const Float_t kA = 0.8303; + const Float_t kB = -0.00392; + const Float_t kC = 0.472 * 0.472; + const Float_t kD = 2.19; - Float_t pr = kA * (kB + TMath::Exp(-x*x / (2. * kC2))); + Float_t pr = kA * (kB + TMath::Exp(-TMath::Power(x*x,kD) / (2.*kC))); return (pr); diff --git a/TRD/AliTRDclusterizerV1.h b/TRD/AliTRDclusterizerV1.h index bd46599df96..4bc178734f8 100644 --- a/TRD/AliTRDclusterizerV1.h +++ b/TRD/AliTRDclusterizerV1.h @@ -28,28 +28,25 @@ class AliTRDclusterizerV1 : public AliTRDclusterizer { virtual Bool_t MakeClusters(); virtual Bool_t ReadDigits(); - virtual void SetClusMaxThresh(Float_t thresh) { fClusMaxThresh = thresh; }; - virtual void SetClusSigThresh(Float_t thresh) { fClusSigThresh = thresh; }; - virtual void SetClusMethod(Int_t meth) { fClusMethod = meth; }; + virtual void SetClusMaxThresh(Int_t thresh) { fClusMaxThresh = thresh; }; + virtual void SetClusSigThresh(Int_t thresh) { fClusSigThresh = thresh; }; - virtual Float_t GetClusMaxThresh() const { return fClusMaxThresh; }; - virtual Float_t GetClusSigThresh() const { return fClusSigThresh; }; - virtual Int_t GetClusMethod() const { return fClusMethod; }; + virtual Int_t GetClusMaxThresh() const { return fClusMaxThresh; }; + virtual Int_t GetClusSigThresh() const { return fClusSigThresh; }; protected: AliTRDdigitsManager *fDigitsManager; //! TRD digits manager - Float_t fClusMaxThresh; // Threshold value for cluster maximum - Float_t fClusSigThresh; // Threshold value for cluster signal - Int_t fClusMethod; // Clustering method + Int_t fClusMaxThresh; // Threshold value for cluster maximum + Int_t fClusSigThresh; // Threshold value for cluster signal private: virtual Float_t Unfold(Float_t eps, Float_t *padSignal); virtual Float_t PadResponse(Float_t x); - ClassDef(AliTRDclusterizerV1,1) // TRD-Cluster manager, slow simulator + ClassDef(AliTRDclusterizerV1,2) // TRD-Cluster finder, slow simulator }; -- 2.43.0