#include <Riostream.h>
#include <TMath.h>
#include <TTree.h>
#include <TObjArray.h>
#include <TClonesArray.h>
#include <TFile.h>
#include <TBranch.h>
#include <TNtuple.h>
#include <TParticle.h>
#include <TGeoMatrix.h>
#include "AliGeomManager.h"
#include "AliPMDcluster.h"
#include "AliPMDclupid.h"
#include "AliPMDrecpoint1.h"
#include "AliPMDrecdata.h"
#include "AliPMDrechit.h"
#include "AliPMDUtility.h"
#include "AliPMDDiscriminator.h"
#include "AliPMDEmpDiscriminator.h"
#include "AliPMDtracker.h"
#include "AliESDPmdTrack.h"
#include "AliESDEvent.h"
#include "AliLog.h"
ClassImp(AliPMDtracker)
AliPMDtracker::AliPMDtracker():
fTreeR(0),
fRecpoints(new TClonesArray("AliPMDrecpoint1", 10)),
fRechits(new TClonesArray("AliPMDrechit", 10)),
fPMDcontin(new TObjArray()),
fPMDcontout(new TObjArray()),
fPMDutil(new AliPMDUtility()),
fPMDrecpoint(0),
fPMDclin(0),
fPMDclout(0),
fXvertex(0.),
fYvertex(0.),
fZvertex(0.),
fSigmaX(0.),
fSigmaY(0.),
fSigmaZ(0.)
{
}
AliPMDtracker:: AliPMDtracker(const AliPMDtracker & ):
TObject(),
fTreeR(0),
fRecpoints(NULL),
fRechits(NULL),
fPMDcontin(NULL),
fPMDcontout(NULL),
fPMDutil(NULL),
fPMDrecpoint(0),
fPMDclin(0),
fPMDclout(0),
fXvertex(0.),
fYvertex(0.),
fZvertex(0.),
fSigmaX(0.),
fSigmaY(0.),
fSigmaZ(0.)
{
AliError("Copy constructor not allowed");
}
AliPMDtracker& AliPMDtracker::operator=(const AliPMDtracker & )
{
AliError("Assignment operator not allowed");
return *this;
}
AliPMDtracker::~AliPMDtracker()
{
if (fRecpoints)
{
fRecpoints->Clear();
}
if (fRechits)
{
fRechits->Clear();
}
if (fPMDcontin)
{
fPMDcontin->Delete();
delete fPMDcontin;
fPMDcontin=0;
}
if (fPMDcontout)
{
fPMDcontout->Delete();
delete fPMDcontout;
fPMDcontout=0;
}
delete fPMDutil;
}
void AliPMDtracker::LoadClusters(TTree *treein)
{
fTreeR = treein;
}
void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
{
Int_t idet = 0;
Int_t ismn = 0;
Int_t trackno = 1, trackpid = 0;
Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
Int_t *irow;
Int_t *icol;
Int_t *itra;
Int_t *ipid;
Float_t *cadc;
AliPMDrechit *rechit = 0x0;
TBranch *branch = fTreeR->GetBranch("PMDRecpoint");
if (!branch)
{
AliError("PMDRecpoint branch not found");
return;
}
branch->SetAddress(&fRecpoints);
TBranch *branch1 = fTreeR->GetBranch("PMDRechit");
if (!branch1)
{
AliError("PMDRechit branch not found");
return;
}
branch1->SetAddress(&fRechits);
Int_t ncrhit = 0;
Int_t nmodules = (Int_t) branch->GetEntries();
AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
for (Int_t imodule = 0; imodule < nmodules; imodule++)
{
branch->GetEntry(imodule);
Int_t nentries = fRecpoints->GetLast();
AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
,nentries));
for(Int_t ient = 0; ient < nentries+1; ient++)
{
fPMDrecpoint = (AliPMDrecpoint1*)fRecpoints->UncheckedAt(ient);
idet = fPMDrecpoint->GetDetector();
ismn = fPMDrecpoint->GetSMNumber();
clusdata[0] = fPMDrecpoint->GetClusX();
clusdata[1] = fPMDrecpoint->GetClusY();
clusdata[2] = fPMDrecpoint->GetClusADC();
clusdata[3] = fPMDrecpoint->GetClusCells();
clusdata[4] = fPMDrecpoint->GetClusSigmaX();
clusdata[5] = fPMDrecpoint->GetClusSigmaY();
if (clusdata[4] >= 0. && clusdata[5] >= 0.)
{
branch1->GetEntry(ncrhit);
Int_t nenbr1 = fRechits->GetLast() + 1;
irow = new Int_t[nenbr1];
icol = new Int_t[nenbr1];
itra = new Int_t[nenbr1];
ipid = new Int_t[nenbr1];
cadc = new Float_t[nenbr1];
for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
{
irow[ient1] = -99;
icol[ient1] = -99;
itra[ient1] = -99;
ipid[ient1] = -99;
cadc[ient1] = 0.;
}
for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
{
rechit = (AliPMDrechit*)fRechits->UncheckedAt(ient1);
itra[ient1] = rechit->GetCellTrack();
ipid[ient1] = rechit->GetCellPid();
cadc[ient1] = rechit->GetCellAdc();
}
if (idet == 0)
{
AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
trackno, trackpid);
}
else if (idet == 1)
{
trackno = itra[0];
trackpid = ipid[0];
}
delete [] irow;
delete [] icol;
delete [] itra;
delete [] ipid;
delete [] cadc;
fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
fPMDcontin->Add(fPMDclin);
ncrhit++;
}
}
}
AliPMDEmpDiscriminator pmddiscriminator;
pmddiscriminator.Discrimination(fPMDcontin,fPMDcontout);
Double_t sectr[4][3] = { {0.,0.,0.},{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
TString snsector="PMD/Sector";
TString symname;
TGeoHMatrix gpmdor;
for(Int_t isector=1; isector<=4; isector++)
{
symname = snsector;
symname += isector;
TGeoHMatrix *gpmdal = AliGeomManager::GetMatrix(symname);
Double_t *tral = gpmdal->GetTranslation();
AliGeomManager::GetOrigGlobalMatrix(symname, gpmdor);
Double_t *tror = gpmdor.GetTranslation();
for(Int_t ixyz=0; ixyz<3; ixyz++)
{
sectr[isector-1][ixyz] = tral[ixyz] - tror[ixyz];
}
}
const Float_t kzpos = 361.5;
Int_t ix = -1, iy = -1;
Int_t det = 0, smn = 0, trno = 1, trpid = 0, mstat = 0;
Float_t xpos = 0., ypos = 0.;
Float_t adc = 0., ncell = 0., radx = 0., rady = 0.;
Float_t xglobal = 0., yglobal = 0., zglobal = 0;
Float_t pid = 0.;
fPMDutil->ApplyAlignment(sectr);
Int_t nentries2 = fPMDcontout->GetEntries();
AliDebug(1,Form("Number of clusters coming after discrimination = %d"
,nentries2));
for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
{
fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
det = fPMDclout->GetDetector();
smn = fPMDclout->GetSMN();
trno = fPMDclout->GetClusTrackNo();
trpid = fPMDclout->GetClusTrackPid();
mstat = fPMDclout->GetClusMatching();
xpos = fPMDclout->GetClusX();
ypos = fPMDclout->GetClusY();
adc = fPMDclout->GetClusADC();
ncell = fPMDclout->GetClusCells();
radx = fPMDclout->GetClusSigmaX();
if ((radx > 999. && radx < 1000.) && ncell == 1)
{
if (smn < 12)
{
ix = (Int_t) (ypos +0.5);
iy = (Int_t) xpos;
}
else if (smn > 12 && smn < 24)
{
ix = (Int_t) xpos;
iy = (Int_t) (ypos +0.5);
}
rady = (Float_t) (ix*100 + iy);
}
else
{
rady = fPMDclout->GetClusSigmaY();
}
pid = fPMDclout->GetClusPID();
if (det == 0)
{
zglobal = kzpos + 1.65;
}
else if (det == 1)
{
zglobal = kzpos - 1.65;
}
fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal,zglobal);
AliESDPmdTrack *esdpmdtr = new AliESDPmdTrack();
esdpmdtr->SetDetector(det);
esdpmdtr->SetSmn(smn);
esdpmdtr->SetClusterTrackNo(trno);
esdpmdtr->SetClusterTrackPid(trpid);
esdpmdtr->SetClusterMatching(mstat);
esdpmdtr->SetClusterX(xglobal);
esdpmdtr->SetClusterY(yglobal);
esdpmdtr->SetClusterZ(zglobal);
esdpmdtr->SetClusterADC(adc);
esdpmdtr->SetClusterCells(ncell);
esdpmdtr->SetClusterPID(pid);
esdpmdtr->SetClusterSigmaX(radx);
esdpmdtr->SetClusterSigmaY(rady);
event->AddPmdTrack(esdpmdtr);
delete esdpmdtr;
}
fPMDcontin->Delete();
fPMDcontout->Delete();
}
void AliPMDtracker::AssignTrPidToCluster(Int_t nentry, Int_t *itra,
Int_t *ipid, Float_t *cadc,
Int_t &trackno, Int_t &trackpid)
{
Int_t *phentry = new Int_t [nentry];
Int_t *hadentry = new Int_t [nentry];
Int_t *trpid = 0x0;
Int_t *sortcoord = 0x0;
Float_t *trenergy = 0x0;
Int_t ngtrack = 0;
Int_t nhtrack = 0;
for (Int_t i = 0; i < nentry; i++)
{
phentry[i] = -1;
hadentry[i] = -1;
if (ipid[i] == 22)
{
phentry[ngtrack] = i;
ngtrack++;
}
else if (ipid[i] != 22)
{
hadentry[nhtrack] = i;
nhtrack++;
}
}
Int_t nghadtrack = ngtrack + nhtrack;
if (ngtrack == 0)
{
trackpid = 8;
trackno = -1;
}
else if (ngtrack >= 1)
{
trenergy = new Float_t [nghadtrack];
trpid = new Int_t [nghadtrack];
sortcoord = new Int_t [2*nghadtrack];
for (Int_t i = 0; i < ngtrack; i++)
{
trenergy[i] = 0.;
trpid[i] = -1;
for (Int_t j = 0; j < nentry; j++)
{
if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
{
trenergy[i] += cadc[j];
trpid[i] = 22;
}
}
}
for (Int_t i = ngtrack; i < nghadtrack; i++)
{
trenergy[i] = 0.;
trpid[i] = -1;
for (Int_t j = 0; j < nentry; j++)
{
if (ipid[j] != 22 && itra[j] == itra[hadentry[i-ngtrack]])
{
trenergy[i] += cadc[j];
trpid[i] = ipid[j];
}
}
}
Bool_t jsort = true;
TMath::Sort(nghadtrack,trenergy,sortcoord,jsort);
Int_t gtr = sortcoord[0];
if (trpid[gtr] == 22)
{
trackpid = 22;
trackno = itra[phentry[gtr]];
}
else
{
trackpid = 8;
trackno = -1;
}
delete [] trenergy;
delete [] trpid;
delete [] sortcoord;
}
delete [] phentry;
delete [] hadentry;
}
void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
{
fXvertex = vtx[0];
fYvertex = vtx[1];
fZvertex = vtx[2];
fSigmaX = evtx[0];
fSigmaY = evtx[1];
fSigmaZ = evtx[2];
}
void AliPMDtracker::ResetClusters()
{
if (fRecpoints) fRecpoints->Clear();
}