#include <TObjArray.h>
#include <TClonesArray.h>
+#include "AliLog.h"
#include "AliRunLoader.h"
#include "AliLoader.h"
#include "AliRawReader.h"
#include "AliPMDrecpoint1.h"
#include "AliPMDRawStream.h"
+
+
ClassImp(AliPMDClusterFinder)
+AliPMDClusterFinder::AliPMDClusterFinder():
+ fRunLoader(0),
+ fPMDLoader(0),
+ fTreeD(0),
+ fTreeR(0),
+ fDigits(new TClonesArray("AliPMDdigit", 1000)),
+ fRecpoints(new TClonesArray("AliPMDrecpoint1", 1000)),
+ fNpoint(0),
+ fEcut(0.)
+{
+//
+// Constructor
+//
+}
+// ------------------------------------------------------------------------- //
AliPMDClusterFinder::AliPMDClusterFinder(AliRunLoader* runLoader):
fRunLoader(runLoader),
fPMDLoader(runLoader->GetLoader("PMDLoader")),
fDigits(new TClonesArray("AliPMDdigit", 1000)),
fRecpoints(new TClonesArray("AliPMDrecpoint1", 1000)),
fNpoint(0),
- fDebug(0),
fEcut(0.)
{
//
TObjArray *pmdcont = new TObjArray();
AliPMDClustering *pmdclust = new AliPMDClustering();
- pmdclust->SetDebug(fDebug);
+
pmdclust->SetEdepCut(fEcut);
fRunLoader->GetEvent(ievt);
- //cout << " ***** Beginning::Digits2RecPoints *****" << endl;
+
fTreeD = fPMDLoader->TreeD();
if (fTreeD == 0x0)
{
- cout << " Can not get TreeD" << endl;
+ AliFatal("AliPMDClusterFinder: Can not get TreeD");
+
}
AliPMDdigit *pmddigit;
TBranch *branch = fTreeD->GetBranch("PMDDigit");
pmdclust->DoClust(idet,ismn,fCellADC,pmdcont);
Int_t nentries1 = pmdcont->GetEntries();
-// cout << " nentries1 = " << nentries1 << endl;
+
+ AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
+
for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
{
AliPMDcluster *pmdcl = (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
delete pmdclust;
delete pmdcont;
- // cout << " ***** End::Digits2RecPoints *****" << endl;
+}
+// ------------------------------------------------------------------------- //
+
+void AliPMDClusterFinder::Digits2RecPoints(AliRawReader *rawReader,
+ TTree *clustersTree)
+{
+ // Converts RAW data to recpoints after running clustering
+ // algorithm on CPV and PREshower plane
+ //
+
+ Float_t clusdata[5];
+
+ TObjArray *pmdcont = new TObjArray();
+ AliPMDClustering *pmdclust = new AliPMDClustering();
+
+ pmdclust->SetEdepCut(fEcut);
+
+ ResetRecpoint();
+
+ Int_t bufsize = 16000;
+ clustersTree->Branch("PMDRecpoint", &fRecpoints, bufsize);
+
+ const Int_t kDDL = 6;
+ const Int_t kRow = 48;
+ const Int_t kCol = 96;
+
+ Int_t idet = 0;
+ Int_t iSMN = 0;
+
+ for (Int_t indexDDL = 0; indexDDL < kDDL; indexDDL++)
+ {
+ if (indexDDL < 4)
+ {
+ iSMN = 6;
+ }
+ else if (indexDDL >= 4)
+ {
+ iSMN = 12;
+ }
+ Int_t ***precpvADC;
+ precpvADC = new int **[iSMN];
+ for (Int_t i=0; i<iSMN; i++) precpvADC[i] = new int *[kRow];
+ for (Int_t i=0; i<iSMN;i++)
+ {
+ for (Int_t j=0; j<kRow; j++) precpvADC[i][j] = new int [kCol];
+ }
+ for (Int_t i = 0; i < iSMN; i++)
+ {
+ for (Int_t j = 0; j < kRow; j++)
+ {
+ for (Int_t k = 0; k < kCol; k++)
+ {
+ precpvADC[i][j][k] = 0;
+ }
+ }
+ }
+ ResetCellADC();
+ rawReader->Reset();
+ AliPMDRawStream pmdinput(rawReader);
+ rawReader->Select(12, indexDDL, indexDDL);
+ while(pmdinput.Next())
+ {
+ Int_t det = pmdinput.GetDetector();
+ Int_t smn = pmdinput.GetSMN();
+ //Int_t mcm = pmdinput.GetMCM();
+ //Int_t chno = pmdinput.GetChannel();
+ Int_t row = pmdinput.GetRow();
+ Int_t col = pmdinput.GetColumn();
+ Int_t sig = pmdinput.GetSignal();
+
+ Int_t indexsmn = 0;
+
+ if (indexDDL < 4)
+ {
+ if (det != 0)
+ AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
+ indexDDL, det));
+ indexsmn = smn - indexDDL * 6;
+ }
+ else if (indexDDL == 4)
+ {
+ if (det != 1)
+ AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
+ indexDDL, det));
+ if (smn < 6)
+ {
+ indexsmn = smn;
+ }
+ else if (smn >= 12 && smn < 18)
+ {
+ indexsmn = smn - 6;
+ }
+ }
+ else if (indexDDL == 5)
+ {
+ if (det != 1)
+ AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
+ indexDDL, det));
+ if (smn >= 6 && smn < 12)
+ {
+ indexsmn = smn - 6;
+ }
+ else if (smn >= 18 && smn < 24)
+ {
+ indexsmn = smn - 12;
+ }
+ }
+ precpvADC[indexsmn][row][col] = sig;
+ } // while loop
+
+ Int_t ismn = 0;
+ for (Int_t indexsmn = 0; indexsmn < iSMN; indexsmn++)
+ {
+ ResetCellADC();
+ for (Int_t irow = 0; irow < kRow; irow++)
+ {
+ for (Int_t icol = 0; icol < kCol; icol++)
+ {
+ fCellADC[irow][icol] =
+ (Double_t) precpvADC[indexsmn][irow][icol];
+ } // row
+ } // col
+ if (indexDDL < 4)
+ {
+ ismn = indexsmn + indexDDL * 6;
+ idet = 0;
+ }
+ else if (indexDDL == 4)
+ {
+ if (indexsmn < 6)
+ {
+ ismn = indexsmn;
+ }
+ else if (indexsmn >= 6 && indexsmn < 12)
+ {
+ ismn = indexsmn + 6;
+ }
+ idet = 1;
+ }
+ else if (indexDDL == 5)
+ {
+ if (indexsmn < 6)
+ {
+ ismn = indexsmn + 6;
+ }
+ else if (indexsmn >= 6 && indexsmn < 12)
+ {
+ ismn = indexsmn + 12;
+ }
+ idet = 1;
+ }
+
+
+ pmdclust->DoClust(idet,ismn,fCellADC,pmdcont);
+ Int_t nentries1 = pmdcont->GetEntries();
+
+ AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
+
+ for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
+ {
+ AliPMDcluster *pmdcl =
+ (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
+ idet = pmdcl->GetDetector();
+ ismn = pmdcl->GetSMN();
+ clusdata[0] = pmdcl->GetClusX();
+ clusdata[1] = pmdcl->GetClusY();
+ clusdata[2] = pmdcl->GetClusADC();
+ clusdata[3] = pmdcl->GetClusCells();
+ clusdata[4] = pmdcl->GetClusRadius();
+
+ AddRecPoint(idet,ismn,clusdata);
+ }
+ pmdcont->Clear();
+
+ //fTreeR->Fill();
+ clustersTree->Fill();
+ ResetRecpoint();
+
+
+ } // smn
+
+ for (Int_t i=0; i<iSMN; i++)
+ {
+ for (Int_t j=0; j<kRow; j++) delete [] precpvADC[i][j];
+ }
+ for (Int_t i=0; i<iSMN; i++) delete [] precpvADC[i];
+ delete precpvADC;
+ } // DDL Loop
+
+ ResetCellADC();
+
+ // delete the pointers
+ delete pmdclust;
+ delete pmdcont;
+
}
// ------------------------------------------------------------------------- //
TObjArray *pmdcont = new TObjArray();
AliPMDClustering *pmdclust = new AliPMDClustering();
- pmdclust->SetDebug(fDebug);
+
pmdclust->SetEdepCut(fEcut);
fRunLoader->GetEvent(ievt);
if (indexDDL < 4)
{
if (det != 0)
- printf(" *** DDL %d and Detector NUMBER %d NOT MATCHING *** ",
- indexDDL, det);
+ AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
+ indexDDL, det));
indexsmn = smn - indexDDL * 6;
}
else if (indexDDL == 4)
{
if (det != 1)
- printf(" *** DDL %d and Detector NUMBER %d NOT MATCHING *** ",
- indexDDL, det);
+ AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
+ indexDDL, det));
if (smn < 6)
{
indexsmn = smn;
else if (indexDDL == 5)
{
if (det != 1)
- printf(" *** DDL %d and Detector NUMBER %d NOT MATCHING *** ",
- indexDDL, det);
+ AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
+ indexDDL, det));
if (smn >= 6 && smn < 12)
{
indexsmn = smn - 6;
pmdclust->DoClust(idet,ismn,fCellADC,pmdcont);
Int_t nentries1 = pmdcont->GetEntries();
- // cout << " nentries1 = " << nentries1 << endl;
+
+ AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
+
for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
{
AliPMDcluster *pmdcl =
delete pmdclust;
delete pmdcont;
- // cout << " ***** End::Digits2RecPoints :: Raw *****" << endl;
}
// ------------------------------------------------------------------------- //
void AliPMDClusterFinder::SetCellEdepCut(Float_t ecut)
fEcut = ecut;
}
// ------------------------------------------------------------------------- //
-void AliPMDClusterFinder::SetDebug(Int_t idebug)
-{
- fDebug = idebug;
-}
-// ------------------------------------------------------------------------- //
void AliPMDClusterFinder::AddRecPoint(Int_t idet,Int_t ismn,Float_t *clusdata)
{
// Add Reconstructed points
#include "Riostream.h"
#include <TNtuple.h>
#include <TObjArray.h>
+#include <stdio.h>
+
#include "AliPMDcluster.h"
#include "AliPMDClustering.h"
-#include <stdio.h>
+#include "AliLog.h"
ClassImp(AliPMDClustering)
const Double_t AliPMDClustering::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
AliPMDClustering::AliPMDClustering():
- fDebug(0),
fCutoff(0.0)
{
for(int i = 0; i < kNDIMX; i++)
if (fEdepCell[i1][i2] > cutoff ) nmx1 = nmx1 + 1;
}
// nmx1 --- number of cells having ener dep >= cutoff
- if (fDebug == 1)
- {
- cout << " nmx1 " << nmx1 << endl;
- }
+
+ AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
// if (nmx1 == 0 | nmx1 == -1) return;
if (nmx1 == 0) nmx1 = 1;
ave=ave/nmx1;
- if (fDebug == 1)
- {
- cout <<"kNMX " << kNMX << " nmx1 " << nmx1<< " ave "<<ave<<
- " cutoff " << cutoff << endl;
- }
+ AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
+ kNMX,ave));
+
incr = CrClust(ave, cutoff, nmx1);
RefClust(incr);
- if (fDebug == 1)
- {
- cout << "fClno " << fClno << endl;
- }
+ AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, fClno));
+
for(i1=0; i1<=fClno; i1++)
{
Float_t cluXC = (Float_t) fClusters[0][i1];
// ofstream ofl0("cells_loc",ios::out);
// initialize fInfocl[2][kNDIMX][kNDIMY]
- if (fDebug == 1)
- {
- printf(" *** Inside CrClust ** kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f\n",
- kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff);
- }
+ AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
+
for (j=0; j < kNDIMX; j++){
for(k=0; k < kNDIMY; k++){
fInfocl[0][j][k] = 0;
for(i=0; i<incr; i++){
if(fInfcl[0][i] != nsupcl){ nsupcl=nsupcl+1; }
if (nsupcl > 4500) {
- Error("RefClust", "Too many superclusters!");
+ AliWarning("RefClust: Too many superclusters!");
nsupcl = 4500;
break;
}
ncl[nsupcl]=ncl[nsupcl]+1;
}
- if (fDebug == 1)
- {
- cout << " # of cells " <<incr+1 << " # of superclusters " << nsupcl+1
- << endl;
- }
+
+ AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
+ incr+1,nsupcl+1));
+
id=-1;
icl=-1;
for(i=0; i<nsupcl; i++){
// cluster center at the centyer of the cell
// cluster radius = half cell dimension
if (fClno >= 5000) {
- Error("RefClust", "Too many clusters!");
+ AliWarning("RefClust: Too many clusters! more than 5000");
return;
}
fClno = fClno + 1;
id = id + 1;
icl = icl+1;
if (fClno >= 5000) {
- Error("RefClust", "Too many clusters!");
+ AliWarning("RefClust: Too many clusters! more than 5000");
return;
}
fClno = fClno+1;
}
for(j=0; j<=ig; j++){
if (fClno >= 5000) {
- Error("RefClust", "Too many clusters!");
+ AliWarning("RefClust: Too many clusters! more than 5000");
return;
}
fClno = fClno + 1;
cm = 16777213./16777216.;
}
else{
- cout << " wrong initialization " << endl;
+ AliWarning("Wrong initialization");
}
}
else{
fCutoff = decut;
}
// ------------------------------------------------------------------------ //
-void AliPMDClustering::SetDebug(Int_t idebug)
-{
- fDebug = idebug;
-}
AliPMDClustering();
virtual ~AliPMDClustering();
- void DoClust(Int_t idet, Int_t ismn, Double_t celladc[][96],
- TObjArray *pmdcont);
- void Order();
-
- Int_t CrClust(Double_t ave, Double_t cutoff, Int_t nmx1);
- void RefClust(Int_t incr);
- void GaussFit(Int_t ncell, Int_t nclust, Double_t &x,
- Double_t &y, Double_t &z, Double_t &xc,
- Double_t &yc, Double_t &zc, Double_t &rc);
+ void DoClust(Int_t idet, Int_t ismn, Double_t celladc[][96],
+ TObjArray *pmdcont);
+ void Order();
+
+ Int_t CrClust(Double_t ave, Double_t cutoff, Int_t nmx1);
+ void RefClust(Int_t incr);
+ void GaussFit(Int_t ncell, Int_t nclust, Double_t &x,
+ Double_t &y, Double_t &z, Double_t &xc,
+ Double_t &yc, Double_t &zc, Double_t &rc);
Double_t Distance(Double_t x1, Double_t y1,
Double_t x2, Double_t y2);
Double_t Ranmar() const;
- void SetEdepCut(Float_t decut);
- void SetDebug(Int_t idebug);
-
+ void SetEdepCut(Float_t decut);
+
protected:
static const Double_t fgkSqroot3by2; // fgkSqroot3by2 = sqrt(3.)/2.
coord --- x and y coordinates of center of each cell
*/
- Int_t fDebug; // Switch for debug (1:Print, 0:Noprint)
Float_t fCutoff; // Energy(ADC) cutoff per cell before clustering
- ClassDef(AliPMDClustering,3) // Does clustering for PMD
+ ClassDef(AliPMDClustering,5) // Does clustering for PMD
};
#endif