/**************************************************************************
* 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. *
**************************************************************************/
///////////////////////////////////////////////////////////////////////////////
// //
// This class contains the points for the ALICE event display //
// //
//Begin_Html
/*
*/
//End_Html
// //
// //
///////////////////////////////////////////////////////////////////////////////
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "AliRun.h"
#include "AliITSpoints.h"
#include "AliITSdisplay.h"
#include "AliITSMap.h"
#include "AliITSRecPointSDDnew.h"
//const Int_t MAXNipx=1026, MAXNipy=1026;
const Int_t MAXNipx=752, MAXNipy=256;
static AliITSMapA2 *sMap = 0;
static Int_t sModule=0;
ClassImp(AliITSpoints)
//_____________________________________________________________________________
AliITSpoints::AliITSpoints()
{
//
// Default constructor
//
fModuleIndex = -1;
fHitIndex = -1;
fTrackIndex = -1;
fDigitIndex = -1;
fMarker[0] = fMarker[1] = fMarker[2]=0;
fMatrix = 0;
//fConnect=kFALSE;
}
//_____________________________________________________________________________
AliITSpoints::AliITSpoints(Int_t npoints)
:AliPoints(npoints)
{
//
// Standard constructor
//
fModuleIndex = -1;
fHitIndex = -1;
fTrackIndex = -1;
fDigitIndex = -1;
fMarker[0] = fMarker[1] = fMarker[2]=0;
fMatrix = 0;
//fConnect=kFALSE;
}
//_____________________________________________________________________________
AliITSpoints::~AliITSpoints()
{
//
// Default destructor
//
fHitIndex = 0;
fTrackIndex = 0;
fDigitIndex = 0;
fModuleIndex = 0;
for (Int_t i=0;i<3;i++){
if (fMarker[i]) delete fMarker[i];
}
fMatrix = 0;
}
//_____________________________________________________________________________
void AliITSpoints::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{
//
//*-*-*-*-*-*-*-*-*-*Execute action corresponding to one event*-*-*-*-*-*-*-*
//*-* =========================================
//*-*
//*-* This member function must be implemented to realize the action
//*-* corresponding to the mouse click on the object in the window
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
if (!gPad) return;
gPad->SetCursor(kCross);
TObject *select = gPad->GetSelected();
if(!select) {printf("no select \n"); return;}
if (!select->InheritsFrom("AliITSpoints")) {gPad->SetUniqueID(0); return;}
//erase old position and draw a line at current position
switch (event) {
case kButton1Down:
return;
case kButton1Motion:
AnodeProjection(px,py);
return;
case kButton1Up:
return;
case kButton1Double:
TimeProjection(px,py);
return;
/* why button2 does not work ?
it does not know kButton2Down, kButton2Motion, kButton2Up
case kButton2Down:
printf("Button2Down \n");
return;
case kButton2Motion:
printf("Button2Motion \n");
TimeProjection(px,py);
return;
case kButton2Up:
return;
*/
}
}
//_____________________________________________________________________________
void AliITSpoints::AnodeProjection(Int_t px, Int_t py)
{
//if(!fDrawHist) return;
gPad->SetCursor(kMove);
gPad->GetCanvas()->FeedbackMode(kTRUE);
Float_t xmin = gPad->GetX1();
Float_t xmax = gPad->GetX2();
Float_t ymin = gPad->GetY1();
Float_t ymax = gPad->GetY2();
//printf("xmin,xmax,ymin,ymax %f %f %f %f\n",xmin,xmax,ymin,ymax);
Float_t x,y;
// x,y in the -1,1 scale
x = gPad->AbsPixeltoX(px);
y = gPad->AbsPixeltoY(py);
Float_t scale[3],center[3];
Int_t irep;
gPad->GetView()->FindScope(scale,center,irep);
AliITSdisplay* display=((AliITSdisplay*)(gAlice->Display()));
Int_t anode, timebin;
display->GetPadIxy(anode,timebin,x*scale[0],y*scale[0]);
Int_t amin, amax, tmin, tmax;
display->GetPadIxy(amin,tmin,xmin*scale[0],ymin*scale[0]);
display->GetPadIxy(amax,tmax,xmax*scale[0],ymax*scale[0]);
if (xmin < 0 && xmax > 0) {
if (0-xmin > xmax) amax-=374;
else {amin+=374; tmin=tmax;}
tmax=256;
}
Int_t tminold=tmin;
Int_t tmaxold=tmax;
tmin=TMath::Min(tminold,tmaxold);
tmax=TMath::Max(tminold,tmaxold);
Int_t nbiny = tmax-tmin;
//printf("amin amax tmin tmax anode, timebin %d %d %d %d %d %d\n",amin,amax,tmin,tmax,anode,timebin);
//create or set the new canvas c2
TVirtualPad *padsav = gPad;
//printf("AnodeProj: gPad %p\n",gPad);
TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
if(c2) delete c2->GetPrimitive("Anode Projection");
if(c2) delete c2->GetPrimitive("Time Projection");
else c2=new TCanvas("c2");
c2->cd();
c2->SetFillColor(0);
char title[30];
sprintf(title,"anode=%d",anode);
TH1F *h1=new TH1F("Anode Projection","",nbiny,(float)tmin,(float)tmax);
h1->SetTitle(title);
h1->SetTitleSize(2.);
h1->SetTitleOffset(2.);
h1->SetLabelSize(0.05);
h1->SetYTitle(title);
//h1->SetStats(0);
AliITSMapA2 *sMap=display->GetMap();
for (int j=tmin;jGetSignal(anode-1,j-1));
h1->Fill((float)j,(float)sMap->GetSignal(anode-1,j-1));
}
h1->Smooth(1);
h1->Fit("gaus","ql");
h1->Draw();
c2->Update();
padsav->cd();
//printf("AnodeProj: gPad padsav %p %p\n",gPad,padsav);
}
//_____________________________________________________________________________
void AliITSpoints::TimeProjection(Int_t px, Int_t py)
{
//if(!fDrawHist) return;
//printf("TimeProjection \n");
gPad->SetCursor(kRightSide);
gPad->GetCanvas()->FeedbackMode(kTRUE);
Float_t xmin = gPad->GetX1();
Float_t xmax = gPad->GetX2();
Float_t ymin = gPad->GetY1();
Float_t ymax = gPad->GetY2();
//printf("xmin,xmax,ymin,ymax %f %f %f %f\n",xmin,xmax,ymin,ymax);
Float_t x,y;
// x,y in the -1,1 scale
x = gPad->AbsPixeltoX(px);
y = gPad->AbsPixeltoY(py);
Float_t scale[3],center[3];
Int_t irep;
gPad->GetView()->FindScope(scale,center,irep);
AliITSdisplay* display=((AliITSdisplay*)(gAlice->Display()));
Int_t anode, timebin;
display->GetPadIxy(anode,timebin,x*scale[0],y*scale[2]);
Int_t amin, amax, tmin, tmax;
display->GetPadIxy(amin,tmin,xmin*scale[0],ymin*scale[2]);
display->GetPadIxy(amax,tmax,xmax*scale[0],ymax*scale[2]);
if (xmin < 0 && xmax > 0) {
if (0-xmin > xmax) amax-=374;
else {amin+=374;tmin=tmax;}
tmax=256;
}
Int_t aminold=amin;
Int_t amaxold=amax;
amin=TMath::Min(aminold,amaxold);
amax=TMath::Max(aminold,amaxold);
Int_t nbinx = amax-amin;
//printf("amin amax tmin tmax anode, timebin %d %d %d %d %d %d\n",amin,amax,tmin,tmax,anode,timebin);
//create or set the new canvas c2
TVirtualPad *padsav = gPad;
TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
if(c2) delete c2->GetPrimitive("Anode Projection");
if(c2) delete c2->GetPrimitive("Time Projection");
else c2=new TCanvas("c2");
c2->cd();
c2->SetFillColor(0);
char title[30];
//printf("title: anode, timebin %d %d\n",anode,timebin);
sprintf(title,"time sample=%d",timebin);
TH1F *h2=new TH1F("Time Projection","",nbinx,(float)amin,(float)amax);
h2->SetTitle(title);
h2->SetTitleSize(2.);
h2->SetTitleOffset(2.);
h2->SetLabelSize(0.05);
h2->SetYTitle(title);
//h2->SetStats(0);
AliITSMapA2 *sMap=display->GetMap();
for (int i=amin;iGetSignal(i-1,timebin-1));
h2->Fill((float)i,(float)sMap->GetSignal(i-1,timebin-1));
}
h2->Smooth(1);
h2->Fit("gaus","ql");
h2->Draw();
c2->Update();
padsav->cd();
}
//_____________________________________________________________________________
void AliITSpoints::DisplayModule()
{
/*
if (!((AliITSDisplay*)(gAlice->Display()))->GetDrawTracksOpt()) return;
// ????
((AliITSDisplay*)(gAlice->Display()))->SetModuleNumber(fModuleIndex);
Int_t event=((AliITSDisplay*)(gAlice->Display()))->GetEvent();
((AliITSDisplay*)(gAlice->Display()))->SetEvent(event);
((AliITSDisplay*)(gAlice->Display()))->SetRange(4,4);
//((AliITSDisplay*)(gAlice->Display()))->DrawModule(fModuleIndex);
*/
}
//_____________________________________________________________________________
const Text_t *AliITSpoints::GetName() const
{
//
// Return name of the Geant3 particle corresponding to this point
//
TParticle *particle = GetParticle();
if (particle) return particle->GetName();
else return IsA()->GetName();
//if (!particle) return "Particle";
}
//_____________________________________________________________________________
Text_t *AliITSpoints::GetObjectInfo(Int_t px, Int_t py)
{
//
// Redefines TObject::GetObjectInfo.
// Displays the info (particle,etc
// corresponding to cursor position px,py
//
if (!gPad) return (char*)"";
static char info[64];
char an[6], tbin[9], track[6];
an="Anode"; tbin="Time bin"; track="Track";
if(strcmp(GetName(),"AliITSpoints ")) {
AliITSdigitSDD *dig=GetDigit();
if (!dig) {sprintf(info,"%s %s %d",GetName(),track,fIndex); return info;}
int anode=dig->fCellX;
int time=dig->fCellY;
sprintf(info,"%s %d %s %d %s %d",an,anode+1,tbin,time+1,track,fIndex);
} else {
if(strcmp(GetName(),"TView ")) {
Float_t x = gPad->AbsPixeltoX(px);
Float_t y = gPad->AbsPixeltoY(py);
sprintf(info,"%s x=%.3g, y=%.3g",GetName(),gPad->PadtoX(x),gPad->PadtoY(y));
} else {
sprintf(info,"%s %s %d",GetName(),track,fIndex);
}
}
return info;
}
//_____________________________________________________________________________
void AliITSpoints::DumpHit()
{
//
// Dump hit corresponding to this point
//
AliITShit *hit = GetHit();
if (hit) hit->Dump();
}
//_____________________________________________________________________________
void AliITSpoints::DumpDigit()
{
//
// Dump digit corresponding to this point
//
AliITSdigitSDD *digit = GetDigit();
if (digit) digit->Dump();
}
//_____________________________________________________________________________
void AliITSpoints::InspectHit()
{
//
// Inspect hit corresponding to this point
//
if (fHitIndex < 0 ) return;
TVirtualPad *padsav = gPad;
AliITShit *hit = GetHit();
if (hit) hit->Inspect();
TVirtualPad *padinspect = (TVirtualPad*)(gROOT->GetListOfCanvases())->FindObject("inspect");
padinspect->cd();
Float_t xmin = gPad->GetX1();
Float_t xmax = gPad->GetX2();
Float_t ymin = gPad->GetY1();
Float_t ymax = gPad->GetY2();
Float_t dy = ymax-ymin;
TPaveText *pad = new TPaveText(xmin, ymin+0.1*dy, xmax, ymin+0.15*dy);
pad->SetBit(kCanDelete);
pad->SetFillColor(42);
pad->Draw();
char ptitle[100];
sprintf(ptitle," %s , fTrack: %d fTrackIndex: %d ",GetName(),fIndex,fTrackIndex);
pad->AddText(ptitle);
padinspect->cd();
padinspect->Update();
if (padsav) padsav->cd();
}
//_____________________________________________________________________________
void AliITSpoints::InspectDigit()
{
//
// Inspect digit corresponding to this point
//
if (fDigitIndex < 0) return;
TVirtualPad *padsav = gPad;
AliITSdigitSDD *digit = GetDigit();
if (digit) digit->Inspect();
TVirtualPad *padinspect = (TVirtualPad*)(gROOT->GetListOfCanvases())->FindObject("inspect");
padinspect->cd();
Float_t xmin = gPad->GetX1();
Float_t xmax = gPad->GetX2();
Float_t ymin = gPad->GetY1();
Float_t ymax = gPad->GetY2();
Float_t dy = ymax-ymin;
TPaveText *pad = new TPaveText(xmin, ymin+0.1*dy, xmax, ymin+0.25*dy);
pad->SetBit(kCanDelete);
pad->SetFillColor(42);
pad->Draw();
char ptitle[11][100];
// sprintf(ptitle[11],"Tracks making this digit");
// pad->AddText(ptitle[11]);
for (int i=0;i<3;i++) {
if (digit->fTracks[i] == 0) continue;
sprintf(ptitle[i],"fTrackIndex: %d Charge: %f",digit->fTracks[i],digit->fTcharges[i]);
pad->AddText(ptitle[i]);
}
padinspect->cd();
padinspect->Update();
if (padsav) padsav->cd();
}
//_____________________________________________________________________________
Int_t AliITSpoints::GetTrackIndex()
{
//
// Dump digit corresponding to this point
//
this->Inspect();
/*
if (fDigitIndex != 0) {
Int_t ncol=this->fMatrix->GetNcols();
for (int i=0;ifMatrix))(0,i),(*(this->fMatrix))(1,i));
}
}
*/
return fTrackIndex;
}
//_____________________________________________________________________________
AliITShit *AliITSpoints::GetHit() const
{
//
// Returns pointer to hit index in AliRun::fParticles
//
AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
gAlice->TreeH()->GetEvent(fTrackIndex);
TClonesArray *ITShits = ITS->Hits();
Int_t nhits = ITShits->GetEntriesFast();
if (fHitIndex < 0 || fHitIndex >= nhits) return 0;
return (AliITShit*)ITShits->UncheckedAt(fHitIndex);
}
//_____________________________________________________________________________
AliITSdigitSDD *AliITSpoints::GetDigit() const
{
//
// Returns pointer to digit index in AliRun::fParticles
//
AliITSdisplay *display=(AliITSdisplay*)gAlice->Display();
Int_t module=display->GetModule();
Int_t layer=display->GetLayer();
Int_t id=(Int_t)((layer-1)/2);
AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
TClonesArray *ITSdigits = ITS->DigitsAddress(id);
gAlice->TreeD()->GetEvent(module);
Int_t ndigits = ITSdigits->GetEntriesFast();
if (fDigitIndex < 0 || fDigitIndex >= ndigits) return 0;
return (AliITSdigitSDD*)ITSdigits->UncheckedAt(fDigitIndex);
/*
// have smth like ??
AliITSgeom *geom = ITS->GetITSgeom;
Int_t lastSPD=geom->GetLastSPD();
Int_t lastSDD=geom->GetLastSDD();
Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
if (nent < lastSDD) {
// take the appropriate action ...using lastSPD, lastSDD ...
// or better introduce an Option_t in the display : All, SPD, SDD, SSD
// or combination
printf("Digitisation wasn't done for all det types !\n");
}
*/
// have smth like
// if (module < lastSPD) return (AliITSdigitSPD*)ITSdigits->UncheckedAt(fDigitIndex)
// else (module < lastSDD) return (AliITSdigitSDD*)ITSdigits->UncheckedAt(fDigitIndex) .... ??
//return 0;
}
//_____________________________________________________________________________
struct Bin {
AliITSdigitSDD *dig;
int idx;
int digidx;
Bin() {dig=0; idx=-1; digidx=-1;}
};
struct PreCluster : public AliITSRecPointSDDnew {
AliITSdigitSDD* summit;
int idx;
int npeaks;
int npoly;
float xpoly[300];
float ypoly[300];
float zpoly[300];
PreCluster() : AliITSRecPointSDDnew() {npeaks=npoly=0;
for (int k=0;k<300;k++) {
xpoly[k]=ypoly[k]=zpoly[k]=0;
}
}
};
//_____________________________________________________________________________
//static void FindCluster(int i, int j, Bin bins[][MAXNipy], PreCluster &c)
static void FindCluster(int i, int j, Bin *bins, PreCluster &c, int thresh)
{
//
// Find clusters
//
//Bin& b=bins[i][j];
Bin& b=bins[i*MAXNipy+j];
Int_t q=b.dig->fSignal - thresh;
printf("i j q %d %d %d\n",i,j,q);
// if (b.idx >= 0 && b.idx != c.idx) {
if (b.idx >= 0 && b.idx > c.idx) {
c.idx=b.idx;
c.npeaks++;
}
if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
// get pad coordinates and prepare the up and down steps
Float_t xl[3], dx[3];
AliITSdisplay *display=(AliITSdisplay*)gAlice->Display();
display->GetPadCxy(i,j,xl,dx);
// calculate center of gravity
c.npoly++;
c.fMulDigits++;
c.fDigitsList[c.npoly-1]=b.digidx;
//printf("npoly c.fDigitsList[c.npoly-1] %d %d \n",c.npoly,c.fDigitsList[c.npoly-1]);
if (c.npoly > 300 ) {
printf("FindCluster - npoly >300, npoly %d \n",c.npoly);
c.npoly=300;
}
c.xpoly[c.npoly-1]=xl[0];
c.ypoly[c.npoly-1]=xl[1];
c.zpoly[c.npoly-1]=xl[2];
c.fX += q*xl[0];
c.fZ += q*xl[2];
c.fY = xl[1];
c.fQ += q;
b.dig = 0; b.idx = c.idx;
/*
// left and right
if (i && bins[i-1][j].dig) FindCluster(i-1,j,bins,c);
if (i < MAXNipx && bins[i+1][j].dig) FindCluster(i+1,j,bins,c);
// up and down
if (j+1 < MAXNipy && bins[i][j+1].dig) FindCluster(i,j+1,bins,c);
if (j && bins[i][j-1].dig) FindCluster(i,j-1,bins,c);
*/
// left and right
if (i && bins[(i-1)*MAXNipy+j].dig) FindCluster(i-1,j,bins,c,thresh);
if (i < MAXNipx && bins[(i+1)*MAXNipy+j].dig) FindCluster(i+1,j,bins,c,thresh);
// up and down
if (j+1 < MAXNipy && bins[i*MAXNipy+(j+1)].dig) FindCluster(i,j+1,bins,c,thresh);
if (j && bins[i*MAXNipy+(j-1)].dig) FindCluster(i,j-1,bins,c,thresh);
}
//_____________________________________________________________________________
void AliITSpoints::GetCenterOfGravity()
{
//
// simple ITS cluster finder from digits -- finds neighbours and
// calculates center of gravity for the cluster
//
const int THRESHOLD=20;
Bin *sBin=new Bin[MAXNipx*MAXNipy];
//else memset(sBin,0,sizeof(Bin)*MAXNipx*MAXNipy);
//printf("sBin %p\n",sBin);
//Bin bins[MAXNipx][MAXNipy];
AliITSdisplay *display=(AliITSdisplay*)gAlice->Display();
Int_t module=display->GetModule(); // or module=fModuleIndex;
Int_t layer=display->GetLayer();
Int_t id=(Int_t)((layer-1)/2);
AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
TClonesArray *ITSdigits = ITS->DigitsAddress(id);
Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
//gAlice->TreeD()->GetEvent(nent-nofmodules+module-1);
gAlice->TreeD()->GetEvent(module);
Int_t ndigits = ITSdigits->GetEntriesFast();
if (fDigitIndex < 0 || fDigitIndex >= ndigits) return;
AliITSdigitSDD *dig;
dig=(AliITSdigitSDD*)ITSdigits->UncheckedAt(fDigitIndex);
Int_t ipx=dig->fCellX;
Int_t ipy=dig->fCellY;
//bins[ipx][ipy].dig=dig;
sBin[ipx*MAXNipy+ipy].dig=dig;
//printf("ipx, ipy, Sbin.dig %d %d %p\n",ipx,ipy,sBin[ipx*MAXNipy+ipy].dig);
int ndig;
int ncls=0;
for (ndig=0;ndigUncheckedAt(ndig);
int i=dig->fCellX; int j=dig->fCellY;
//printf("ndig i j %d %d %d \n",ndig,i,j);
//bins[i][j].dig=dig;
sBin[i*MAXNipy+j].dig=dig;
sBin[i*MAXNipy+j].digidx=ndig;
}
//PreCluster c; c.summit=bins[ipx][ipy].dig; c.idx=ncls;
// FindCluster(ipx, ipy, bins, c);
PreCluster c; c.summit=sBin[ipx*MAXNipy+ipy].dig; c.idx=ncls;
FindCluster(ipx, ipy, sBin, c, 0);
if (c.npeaks>1) {
printf("GetCenterOfGravity -- more than one peak");
}
c.fX /= c.fQ;
c.fZ /= c.fQ;
ncls++;
AliITSpoints *points = 0;
points = new AliITSpoints(1);
points->SetMarkerColor(kYellow);
points->SetMarkerStyle(5);
points->SetMarkerSize(1.);
points->SetPoint(0,c.fX,c.fY,c.fZ);
points->SetParticle(-1);
points->Draw();
TPolyLine3D *pline=0;
/*
pline=new TPolyLine3D(c.npoly);
Int_t np=c.npoly;
TVector *xp=new TVector(c.npoly);
TVector *yp=new TVector(c.npoly);
TVector *zp=new TVector(c.npoly);
for (int i=0;iSetPoint(i,(*xp)(i),(*yp)(i),(*zp)(i));
//printf("np, i, xp, yp, zp %d %d %f %f %f \n",np,i,(*xp)(i),(*yp)(i),(*zp)(i));
}
*/
/*
delete xp;
delete yp;
delete zp;
*/
//if (fConnect) {
pline=new TPolyLine3D(c.npoly,c.xpoly,c.ypoly,c.zpoly);
pline->SetLineColor(kWhite);
pline->Draw();
//}
for (int k=0;kDisplay();
Int_t module=display->GetModule(); // or module=fModuleIndex;
Int_t layer=display->GetLayer();
Int_t id=(Int_t)((layer-1)/2);
AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
TClonesArray *ITSdigits = ITS->DigitsAddress(id);
Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
//gAlice->TreeD()->GetEvent(nent-nofmodules+module-1);
gAlice->TreeD()->GetEvent(module);
Int_t ndigits = ITSdigits->GetEntriesFast();
if (fDigitIndex < 0 || fDigitIndex >= ndigits) return;
AliITSdigitSDD *dig;
dig=(AliITSdigitSDD*)ITSdigits->UncheckedAt(fDigitIndex);
Int_t ipx=dig->fCellX;
Int_t ipy=dig->fCellY;
lBin[ipx*MAXNipy+ipy].dig=dig;
int ndig;
int ncls=0;
const int kMaxNeighbours=4;
for (ndig=0;ndigUncheckedAt(ndig);
int i=dig->fCellX; int j=dig->fCellY;
//printf("ndig i j %d %d %d \n",ndig,i,j);
lBin[i*MAXNipy+j].dig=dig;
lBin[i*MAXNipy+j].digidx=ndig;
}
PreCluster c; c.summit=lBin[ipx*MAXNipy+ipy].dig; c.idx=ncls;
FindCluster(ipx, ipy, lBin, c, 0);
printf("finished FindCluster !!!\n");
Int_t fMul=c.npoly;
printf("fMul %d \n",fMul);
AliITSdigitSDD* digt;
int fIx[fMul];
int fIy[fMul];
int fSortIx[fMul];
int fSortIy[fMul];
int fQ[fMul];
int fQini[fMul];
int fIndLocal[fMul];
int fIndSaddle[fMul];
AliITSdigitSDD* fDig[fMul]; // current list of digits
for (int i=0; iUncheckedAt(idx);
fDig[i]=digt;
fIx[i]= digt->fCellX;
fIy[i]= digt->fCellY;
fSortIx[i]= digt->fCellX;
fSortIy[i]= digt->fCellY;
fQini[i] = digt->fSignal;
}
Int_t idx_sort[fMul];
TMath::Sort(fMul,fSortIy,idx_sort,0);
Int_t tmin=fIy[idx_sort[0]];
Int_t tmax=fIy[idx_sort[fMul-1]];
TMath::Sort(fMul,fSortIx,idx_sort,0);
Int_t amin=fIx[idx_sort[0]];
Int_t amax=fIx[idx_sort[fMul-1]];
//printf("tmin tmax amin amax %d %d %d %d\n",tmin,tmax,amin,amax);
AliITSMapA2 *sMap=display->GetMap();
Int_t nt=tmax-tmin;
Int_t nbinx=amax-amin;
float **source = new float* [nbinx+1];
for(int i=0;iGetSignal(i,j);
Qsmooth[i-amin][j-tmin]=sMap->GetSignal(i,j);
//printf("i j Qsmooth signal map %d %d %f %f\n",i,j,Qsmooth[i-amin][j-tmin],sMap->GetSignal(i,j));
temp[j-tmin]=sMap->GetSignal(i,j);
}
h1->SmoothArray(nt+1,temp,1);
for (int j=tmin;j<=tmax;j++) {
Qsmooth[i-amin][j-tmin]=temp[j-tmin];
//printf("i j Qsmooth temp %d %d %f %f\n",i,j,Qsmooth[i-amin][j-tmin],temp[j-tmin]);
}
}
//
// Find local maxima//
int fNLocal=0;
int fNSaddle=0;
Bool_t IsLocal[fMul];
Bool_t IsSaddle[fMul];
Int_t AssocPeak[fMul];
Int_t nn;
Int_t Qneighb[2*kMaxNeighbours];
Int_t X[2*kMaxNeighbours], Y[2*kMaxNeighbours];
for (int i=0; iamax || Y[j]tmax) continue;
Qneighb[j]=(Int_t)Qsmooth[X[j]-amin][Y[j]-tmin];
}
int nlzero=0;
for (int j=0; j fQ[i] || nlzero>=2) { //???
if (Qneighb[j] > fQ[i] ) {
//printf("i j %d %d local kFALSE Qneighb fQ %d %d\n",i,j,Qneighb[j],fQ[i]);
IsLocal[i]=kFALSE;
break;
//
// handle special case of neighbouring pads with equal signal
} else if (Qneighb[j] == fQ[i]) {
if (fNLocal >0) {
for (Int_t k=0; k= tmin && j+2 <= tmax && (fQ[i] == Qsmooth[X[j]-amin][Y[j-2]-tmin]) || (fQ[i] == Qsmooth[X[j]-amin][Y[j+2]-tmin]))
{
IsLocal[i]=kFALSE;
}
} // loop over local maxima
} // are there are already local maxima
} // IsLocal
} // loop over neighb
// find saddle points
int nzero=0; int nlower=0;
for (int j=0; jamax || Y[j]tmax) continue;
Qneighb[j]=(int)source[X[j]-amin][Y[j]-tmin];
if (!Qneighb[j] || !fQ[i] || !fQini[i]) {nzero++;continue;}
//if (!fQ[i] || Qneighb[j]<4 || !fQini[i]) {nzero++;continue;}
if (Qneighb[j] < fQ[i] ) {
//printf("fIx fIy X[j] Y[j] %d %d %d %d saddle kFALSE Qneighb fQ %d %d\n",fIx[i],fIy[i],X[j],Y[j],Qneighb[j],fQ[i]);
nlower++;
IsSaddle[i]=kFALSE;
break;
//
// handle special case of neighbouring pads with equal signal
// this is dangerous for saddle points ! better out !!
} else if (Qneighb[j] == fQ[i]) {
if (fNSaddle >0) {
for (Int_t k=0; k0) {
//printf(" i fNLocal nlzero %d %d %d\n",i,fNLocal,nlzero);
fIndLocal[fNLocal]=i;
fNLocal++;
} else fQ[i]=0;
// Maxima should not be on the edge
//if (IsSaddle[i] && fQ[i]>0 && !nzero ) {
if (IsSaddle[i] && fQini[i]>0) {
//printf(" i fNSaddle nzero %d %d %d\n",i,fNSaddle,nzero);
if (fIy[i]-1 < tmin || fIy[i]+1 > tmax) continue;
if (source[fIx[i]-amin][fIy[i]-1-tmin] && source[fIx[i]-amin][fIy[i]+1-tmin]) {
fIndSaddle[fNSaddle]=i;
fNSaddle++;
}
}
} // loop over all digits
printf("fNLocal %d fNSaddle %d fMul %d \n",fNLocal,fNSaddle,fMul);
// take the highest local maxima and spread the charge, i.e. define the
// response
Int_t idx=TMath::LocMax(fMul,fQ);
int thresh=0;
int qsaddle[fNSaddle];
if (fNSaddle) {
for (int i=0;i2 and should apply only to
// clusters with fNLocal==2
Float_t xc[fNLocal], yc[fNLocal], zc[fNLocal];
if (fNLocal > 1 && fNSaddle) {
for (int i=0; i thresh+3) lBin[fIx[i]*MAXNipy+fIy[i]].dig=fDig[i];
else lBin[fIx[i]*MAXNipy+fIy[i]].dig=0;
printf("fDig %p \n",fDig[i]);
printf("i j source %d %d %f lBin %p \n",fIx[i],fIy[i],source[fIx[i]-amin][fIy[i]-tmin],lBin[fIx[i]*MAXNipy+fIy[i]].dig);
}
for (int i=0; ifSignal);
printf("%p \n",lBin[ipxl*MAXNipy+ipyl].dig);
PreCluster cnew; cnew.summit=lBin[ipxl*MAXNipy+ipyl].dig; cnew.idx=ncls;
printf("q %p %d \n",fDig[fIndLocal[i]],lBin[ipx*MAXNipy+ipy].dig->fSignal);
//FindCluster(ipx, ipy, lBin, cnew, 2*thresh+1); // it's further away
FindCluster(ipxl, ipyl, lBin, cnew, 0);
printf("finish FindCluster i %d\n",i);
cnew.fX /= cnew.fQ;
cnew.fZ /= cnew.fQ;
ncls++;
printf("i cnew.fX cnew.fY cnew.fZ %d %f %f %f\n",i,cnew.fX, cnew.fY, cnew.fZ);
xc[i]=cnew.fX;
yc[i]=cnew.fY;
zc[i]=cnew.fZ;
Int_t fMulnew=cnew.npoly;
printf("fMulnew %d \n",fMulnew);
for (int j=0; jUncheckedAt(idx);
fDig[j]=digt;
//printf("digt %p \n",digt);
fIx[j]= digt->fCellX;
fIy[j]= digt->fCellY;
fSortIx[j]= digt->fCellX;
fSortIy[j]= digt->fCellY;
fQini[j] = digt->fSignal;
// printf("i fIx fIy fQini %d %d %d %d\n",j,fIx[j],fIy[j],fQini[j]);
}
Int_t idx_sort[fMulnew];
TMath::Sort(fMulnew,fSortIy,idx_sort,0);
Int_t tmin=fIy[idx_sort[0]];
Int_t tmax=fIy[idx_sort[fMulnew-1]];
if (ipy-tmin != tmax-ipy) {
printf("non-symetric cluster in time direction\n");
// do something - take Local maxima coord
} else {
// take cog coord
}
printf("ipy-tmin tmax-ipy %d %d\n",ipy-tmin, tmax-ipy);
TMath::Sort(fMulnew,fSortIx,idx_sort,0);
Int_t amin=fIx[idx_sort[0]];
Int_t amax=fIx[idx_sort[fMulnew-1]];
if (ipxl-amin != amax-ipxl) {
printf("non-symetric cluster in anode direction\n");
// do something - take Local maxima coord
} else {
// take cog coord
}
printf("ipxl-amin amax-ipxl %d %d\n",ipxl-amin, amax-ipxl);
printf("inside loop over local maxima: i ipxl ipyl %d %d %d\n",i,ipxl,ipyl);
} // loop over local maxima
} // if fNLocal
if (fNLocal == 1) {
int ipx1=fIx[fIndLocal[0]]; int ipy1=fIy[fIndLocal[0]];
if (ipy1-tmin != tmax-ipy1) {
printf("non-symetric cluster in time direction\n");
// do something - take Local maxima coord or fit
} else {
// take cog coord or fit
}
printf("ipy1-tmin tmax-ipy1 %d %d\n",ipy1-tmin, tmax-ipy1);
if (ipx1-amin != amax-ipx1) {
printf("non-symetric cluster in anode direction\n");
// do something - take Local maxima coord or fit
} else {
// take cog coord or fit
}
printf("ipx1-amin amax-ipx1 %d %d\n",ipx1-amin, amax-ipx1);
}
printf("Here! fNLocal fNSaddle %d %d\n",fNLocal,fNSaddle);
AliITSpoints *points = 0;
for (int i=0; iSetMarkerColor(kYellow);
points->SetMarkerStyle(5);
points->SetMarkerSize(1.);
points->SetPoint(0,xc[i],yc[i],zc[i]);
points->SetParticle(-1);
points->Draw();
}
//points = new AliITSpoints(fNLocal);
for (int i=0; iSetMarkerColor(kGreen);
points->SetMarkerStyle(5);
points->SetMarkerSize(1.5);
points->SetPoint(0,c.xpoly[idx],c.ypoly[idx],c.zpoly[idx]);
points->SetParticle(-1);
points->Draw();
}
//points = new AliITSpoints(fNSaddle);
for (int i=0; iSetMarkerColor(kWhite);
points->SetMarkerStyle(5);
points->SetMarkerSize(1.5);
points->SetPoint(0,c.xpoly[idx],c.ypoly[idx],c.zpoly[idx]);
points->SetParticle(-1);
points->Draw();
}
for (int k=0;k