#include <AliESDVertex.h>
-#include <AliITSVertexer.h>
-#include <AliRunLoader.h>
-#include <AliITSLoader.h>
-#include <AliMultiplicity.h>
-#include <AliITSMultReconstructor.h>
+#include "AliITSgeom.h"
+#include "AliITSVertexer.h"
+#include "AliRunLoader.h"
+#include "AliITSLoader.h"
+#include "AliMultiplicity.h"
+#include "AliITSMultReconstructor.h"
+
+const Float_t AliITSVertexer::fgkPipeRadius = 3.0;
ClassImp(AliITSVertexer)
//////////////////////////////////////////////////////////////////////
//______________________________________________________________________
-AliITSVertexer::AliITSVertexer():AliVertexer() {
+AliITSVertexer::AliITSVertexer():AliVertexer(),
+fLadders(),
+fLadOnLay2(0) {
// Default Constructor
+ SetLaddersOnLayer2();
}
-AliITSVertexer::AliITSVertexer(TString filename):AliVertexer() {
+AliITSVertexer::AliITSVertexer(TString filename):AliVertexer(),
+fLadders(),
+fLadOnLay2(0)
+{
// Standard constructor
AliRunLoader *rl = AliRunLoader::GetRunLoader();
if(!rl){
Fatal("AliITSVertexer","Run Loader not found");
}
+ /*
if(rl->LoadgAlice()){
Fatal("AliITSVertexer","The AliRun object is not available - nothing done");
}
+ */
fCurrentVertex = 0;
SetFirstEvent(0);
SetLastEvent(0);
- rl->LoadHeader();
+ // rl->LoadHeader();
AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
if(!filename.Contains("default"))itsLoader->SetVerticesFileName(filename);
if(!filename.Contains("null"))itsLoader->LoadVertices("recreate");
itsLoader->LoadRecPoints();
- Int_t lst;
+ // Int_t lst;
+ SetLastEvent(rl->GetNumberOfEvents()-1);
+ /*
if(rl->TreeE()){
lst = static_cast<Int_t>(rl->TreeE()->GetEntries());
SetLastEvent(lst-1);
}
+ */
+ SetLaddersOnLayer2();
}
//______________________________________________________________________
return *this;
}
+//______________________________________________________________________
+AliITSVertexer::~AliITSVertexer() {
+ // Destructor
+ if(fLadders) delete [] fLadders;
+}
+
//______________________________________________________________________
void AliITSVertexer::FindMultiplicity(Int_t evnumber){
// Invokes AliITSMultReconstructor to determine the
return;
}
+//______________________________________________________________________
+void AliITSVertexer::SetLaddersOnLayer2(Int_t ladwid){
+ // Calculates the array of ladders on layer 2 to be used with a
+ // given ladder on layer 1
+ fLadOnLay2=ladwid;
+ AliRunLoader *rl =AliRunLoader::GetRunLoader();
+ AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
+ AliITSgeom* geom = itsLoader->GetITSgeom();
+ Int_t ladtot1=geom->GetNladders(1);
+ if(fLadders) delete [] fLadders;
+ fLadders=new UShort_t[ladtot1];
+
+
+ Double_t pos1[3],pos2[3];
+ Int_t mod1=geom->GetModuleIndex(2,1,1);
+ geom->GetTrans(mod1,pos1); // position of the module in the MRS
+ Double_t phi0=TMath::ATan2(pos1[1],pos1[0]);
+ if(phi0<0) phi0+=2*TMath::Pi();
+ Int_t mod2=geom->GetModuleIndex(2,2,1);
+ geom->GetTrans(mod2,pos2);
+ Double_t phi2=TMath::ATan2(pos2[1],pos2[0]);
+ if(phi2<0) phi2+=2*TMath::Pi();
+ Double_t deltaPhi= phi0-phi2; // phi width of a layer2 module
+
+ for(Int_t i= 0; i<ladtot1;i++){
+ Int_t modlad= geom->GetModuleIndex(1,i+1,1);
+ Double_t posmod[3];
+ geom->GetTrans(modlad,posmod);
+ Double_t phimod=TMath::ATan2(posmod[1],posmod[0]);
+ if(phimod<0) phimod+=2*TMath::Pi();
+ Double_t phi1= phimod+deltaPhi*double(fLadOnLay2);
+ if(phi1<0) phi1+=2*TMath::Pi();
+ if(phi1>2*TMath::Pi()) phi1-=2*TMath::Pi();
+ Double_t philad1=phi0-phi1;
+ UShort_t lad1;
+ Double_t ladder1=(philad1)/(deltaPhi) +1.;
+ if(ladder1<1){ladder1=40+ladder1;}
+ lad1=int(ladder1+0.5);
+ fLadders[i]=lad1;
+ }
+}
+
+
//______________________________________________________________________
void AliITSVertexer::WriteCurrentVertex(){
// Write the current AliVertex object to file fOutFile
AliITSVertexer();
// standard constructor
AliITSVertexer(TString filename);
- virtual ~AliITSVertexer(){;}
+ virtual ~AliITSVertexer();
virtual void FindMultiplicity(Int_t evnumber);
virtual void WriteCurrentVertex();
-
+ const Float_t GetPipeRadius()const {return fgkPipeRadius;}
+ virtual void SetLaddersOnLayer2(Int_t ladwid=4);
protected:
// assignment operator (NO assignment allowed)
AliITSVertexer& operator=(const AliITSVertexer& /* vtxr */);
- ClassDef(AliITSVertexer,3);
+ static const Float_t fgkPipeRadius; // beam pipe radius (cm)
+ UShort_t *fLadders; // array with layer1-layer2 ladders correspondances
+ Int_t fLadOnLay2; // (2*fLadOnLay2+1)=number of layer2 ladders
+ // associated to a layer1 ladder
+
+ ClassDef(AliITSVertexer,4);
};
#endif
#include "AliLog.h"
#include "AliRunLoader.h"
#include "AliITSLoader.h"
-#include "AliITSgeom.h"
#include "AliITSDetTypeRec.h"
#include "AliITSRecPoint.h"
/////////////////////////////////////////////////////////////////
fZCutDiamond(0.),
fMaxZCut(0.),
fDCAcut(0.),
-fDiffPhiMax(0.) {
+fDiffPhiMax(0.)
+ {
// Default constructor
SetCoarseDiffPhiCut();
SetCoarseMaxRCut();
SetMaxZCut();
SetDCAcut();
SetDiffPhiMax();
-
}
//______________________________________________________________________
fZCutDiamond(0.),
fMaxZCut(0.),
fDCAcut(0.),
-fDiffPhiMax(0.) {
+fDiffPhiMax(0.)
+{
// Standard constructor
fLines = new TClonesArray("AliStrLine",1000);
SetCoarseDiffPhiCut();
//______________________________________________________________________
AliITSVertexer3D::~AliITSVertexer3D() {
// Destructor
- fVert3D = 0; // this object is not owned by this class
- if(fLines){
+ if(fLines){
fLines->Delete();
delete fLines;
}
}
+//______________________________________________________________________
+void AliITSVertexer3D::ResetVert3D(){
+ //
+ fVert3D.SetXv(0.);
+ fVert3D.SetYv(0.);
+ fVert3D.SetZv(0.);
+ fVert3D.SetDispersion(0.);
+ fVert3D.SetNContributors(0);
+}
//______________________________________________________________________
AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(Int_t evnumber){
// Defines the AliESDVertex for the current event
+ ResetVert3D();
AliDebug(1,Form("FindVertexForCurrentEvent - 3D - PROCESSING EVENT %d",evnumber));
- fVert3D = 0;
if(fLines)fLines->Clear();
Int_t nolines = FindTracklets(evnumber,0);
fCurrentVertex = 0;
if(nolines<2)return fCurrentVertex;
Int_t rc=Prepare3DVertex(0);
- if(rc==0)Find3DVertex();
- if(fVert3D){
+ if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(fLines,0);
+ /* uncomment to debug
+ printf("Vertex found in first iteration:\n");
+ fVert3D.Print();
+ printf("Start second iteration\n");
+ end of debug lines */
+ if(fVert3D.GetNContributors()>0){
if(fLines) fLines->Delete();
nolines = FindTracklets(evnumber,1);
if(nolines>=2){
rc=Prepare3DVertex(1);
- if(rc==0)Find3DVertex();
+ if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(fLines,0);
}
}
-
- if(fVert3D){
+ /* uncomment to debug
+ printf("Vertex found in second iteration:\n");
+ fVert3D.Print();
+ end of debug lines */
+
+ Float_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv());
+ if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){
fCurrentVertex = new AliESDVertex();
fCurrentVertex->SetTitle("vertexer: 3D");
fCurrentVertex->SetName("Vertex");
- fCurrentVertex->SetXv(fVert3D->GetXv());
- fCurrentVertex->SetYv(fVert3D->GetYv());
- fCurrentVertex->SetZv(fVert3D->GetZv());
- fCurrentVertex->SetDispersion(fVert3D->GetDispersion());
- fCurrentVertex->SetNContributors(fVert3D->GetNContributors());
+ fCurrentVertex->SetXv(fVert3D.GetXv());
+ fCurrentVertex->SetYv(fVert3D.GetYv());
+ fCurrentVertex->SetZv(fVert3D.GetZv());
+ fCurrentVertex->SetDispersion(fVert3D.GetDispersion());
+ fCurrentVertex->SetNContributors(fVert3D.GetNContributors());
}
FindMultiplicity(evnumber);
return fCurrentVertex;
Float_t deltaR=fCoarseMaxRCut;
Float_t dZmax=fZCutDiamond;
if(optCuts){
- xbeam=fVert3D->GetXv();
- ybeam=fVert3D->GetYv();
- zvert=fVert3D->GetZv();
+ xbeam=fVert3D.GetXv();
+ ybeam=fVert3D.GetYv();
+ zvert=fVert3D.GetZv();
deltaPhi = fDiffPhiMax;
deltaR=fMaxRCut;
dZmax=fMaxZCut;
Int_t nolines = 0;
// Loop on modules of layer 1
for(Int_t modul1= irstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1
+ UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
branch->GetEvent(modul1);
Int_t nrecp1 = itsRec->GetEntries();
TClonesArray *prpl1 = new TClonesArray("AliITSRecPoint",nrecp1);
geom->LtoG(modul1,lc,gc); // global coordinates
Double_t phi1 = TMath::ATan2(gc[1]-ybeam,gc[0]-xbeam);
if(phi1<0)phi1=2*TMath::Pi()+phi1;
- for(Int_t modul2= irstL2; modul2<=lastL2;modul2++){
- branch->GetEvent(modul2);
- Int_t nrecp2 = itsRec->GetEntries();
- for(Int_t j2=0;j2<nrecp2;j2++){
- recp = (AliITSRecPoint*)itsRec->At(j2);
- lc2[0]=recp->GetDetLocalX();
- lc2[2]=recp->GetDetLocalZ();
- geom->LtoG(modul2,lc2,gc2);
- Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
- if(phi2<0)phi2=2*TMath::Pi()+phi2;
- Double_t diff = TMath::Abs(phi2-phi1);
- if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
- if(diff>deltaPhi)continue;
- AliStrLine line(gc,gc2,kTRUE);
- Double_t cp[3];
- Int_t retcode = line.Cross(&zeta,cp);
- if(retcode<0)continue;
- Double_t dca = line.GetDCA(&zeta);
- if(dca<0.) continue;
- if(dca>deltaR)continue;
- Double_t deltaZ=cp[2]-zvert;
- if(TMath::Abs(deltaZ)>dZmax)continue;
- MakeTracklet(gc,gc2,nolines);
+ for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
+ for(Int_t k=0;k<4;k++){
+ Int_t ladmod=fLadders[ladder-1]+ladl2;
+ if(ladmod>geom->GetNladders(2)) ladmod=ladmod-geom->GetNladders(2);
+ Int_t modul2=geom->GetModuleIndex(2,ladmod,k+1);
+ branch->GetEvent(modul2);
+ Int_t nrecp2 = itsRec->GetEntries();
+ for(Int_t j2=0;j2<nrecp2;j2++){
+ recp = (AliITSRecPoint*)itsRec->At(j2);
+ lc2[0]=recp->GetDetLocalX();
+ lc2[2]=recp->GetDetLocalZ();
+ geom->LtoG(modul2,lc2,gc2);
+ Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
+ if(phi2<0)phi2=2*TMath::Pi()+phi2;
+ Double_t diff = TMath::Abs(phi2-phi1);
+ if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
+ if(diff>deltaPhi)continue;
+ AliStrLine line(gc,gc2,kTRUE);
+ Double_t cp[3];
+ Int_t retcode = line.Cross(&zeta,cp);
+ if(retcode<0)continue;
+ Double_t dca = line.GetDCA(&zeta);
+ if(dca<0.) continue;
+ if(dca>deltaR)continue;
+ Double_t deltaZ=cp[2]-zvert;
+ if(TMath::Abs(deltaZ)>dZmax)continue;
+ MakeTracklet(gc,gc2,nolines);
+ }
+ detTypeRec.ResetRecPoints();
}
- detTypeRec.ResetRecPoints();
}
}
delete prpl1;
}
-//______________________________________________________________________
-void AliITSVertexer3D::Find3DVertex(){
- // Determines the vertex position
- // same algorithm as in AliVertexerTracks::StrLinVertexFinderMinDist(0)
- // adapted to pure AliStrLine objects
-
- Int_t knacc = fLines->GetEntriesFast();
- Double_t (*vectP0)[3]=new Double_t [knacc][3];
- Double_t (*vectP1)[3]=new Double_t [knacc][3];
-
- Double_t sum[3][3];
- Double_t dsum[3]={0,0,0};
- for(Int_t i=0;i<3;i++)
- for(Int_t j=0;j<3;j++)sum[i][j]=0;
- for(Int_t i=0; i<knacc; i++){
- AliStrLine *line1 = (AliStrLine*)fLines->At(i);
-
- Double_t p0[3],cd[3];
- line1->GetP0(p0);
- line1->GetCd(cd);
- Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
- vectP0[i][0]=p0[0];
- vectP0[i][1]=p0[1];
- vectP0[i][2]=p0[2];
- vectP1[i][0]=p1[0];
- vectP1[i][1]=p1[1];
- vectP1[i][2]=p1[2];
-
- Double_t matr[3][3];
- Double_t dknow[3];
- AliVertexerTracks::GetStrLinDerivMatrix(p0,p1,matr,dknow);
-
- for(Int_t iii=0;iii<3;iii++){
- dsum[iii]+=dknow[iii];
- for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
- }
- }
-
- Double_t vett[3][3];
- Double_t det=AliVertexerTracks::GetDeterminant3X3(sum);
- Double_t initPos[3];
- Double_t sigma = 0.;
- for(Int_t i=0;i<3;i++)initPos[i]=0.;
-
- Int_t goodvert=0;
-
- if(det!=0){
- for(Int_t zz=0;zz<3;zz++){
- for(Int_t ww=0;ww<3;ww++){
- for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
- }
- for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
- initPos[zz]=AliVertexerTracks::GetDeterminant3X3(vett)/det;
- }
-
- Double_t rvert=TMath::Sqrt(initPos[0]*initPos[0]+initPos[1]*initPos[1]);
- Double_t zvert=initPos[2];
- if(rvert<fCoarseMaxRCut && TMath::Abs(zvert)<fZCutDiamond){
- for(Int_t i=0; i<knacc; i++){
- Double_t p0[3]={0,0,0},p1[3]={0,0,0};
- for(Int_t ii=0;ii<3;ii++){
- p0[ii]=vectP0[i][ii];
- p1[ii]=vectP1[i][ii];
- }
- sigma+=AliVertexerTracks::GetStrLinMinDist(p0,p1,initPos);
- }
- sigma=TMath::Sqrt(sigma);
- goodvert=1;
- }
- }
- delete vectP0;
- delete vectP1;
- if(!goodvert){
- Warning("Find3DVertex","Finder did not succed");
- for(Int_t i=0;i<3;i++)initPos[i]=0.;
- sigma=999;
- knacc=-1;
- }
- if(fVert3D){
- fVert3D-> SetXYZ(initPos);
- fVert3D->SetDispersion(sigma);
- fVert3D->SetNContributors(knacc);
- }else{
- fVert3D = new AliVertex(initPos,sigma,knacc);
- }
- // fVert3D->Print();
-}
-
-
-
//______________________________________________________________________
Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){
// Finds the 3D vertex information using tracklets
Float_t deltaR=fCoarseMaxRCut;
Float_t dZmax=fZCutDiamond;
if(optCuts){
- xbeam=fVert3D->GetXv();
- ybeam=fVert3D->GetYv();
- zvert=fVert3D->GetZv();
+ xbeam=fVert3D.GetXv();
+ ybeam=fVert3D.GetYv();
+ zvert=fVert3D.GetZv();
deltaR=fMaxRCut;
dZmax=fMaxZCut;
}
Float_t binsizez=(zh-zl)/nbz;
TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh);
-
// cleanup of the TCLonesArray of tracklets (i.e. fakes are removed)
Int_t *validate = new Int_t [fLines->GetEntriesFast()];
for(Int_t i=0; i<fLines->GetEntriesFast();i++)validate[i]=0;
Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY);
if(TMath::Abs(deltaZ)>dZmax)continue;
if(raddist>deltaR)continue;
- validate[i]++;
- validate[j]++;
+ validate[i]=1;
+ validate[j]=1;
h3d->Fill(point[0],point[1],point[2]);
}
}
fLines->Compress();
AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines->GetEntriesFast()));
-
if(fLines->GetEntriesFast()>1){
- Find3DVertex(); // find a first candidate for the primary vertex
+ // find a first candidate for the primary vertex
+ fVert3D=AliVertexerTracks::TrackletVertexFinder(fLines,0);
// make a further selection on tracklets based on this first candidate
- fVert3D->GetXYZ(peak);
- AliDebug(1,Form("FIRSTgv Ma V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2]));
+ fVert3D.GetXYZ(peak);
+ AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2]));
for(Int_t i=0; i<fLines->GetEntriesFast();i++){
AliStrLine *l1 = (AliStrLine*)fLines->At(i);
if(l1->GetDistFromPoint(peak)> fDCAcut)fLines->RemoveAt(i);
}
fLines->Compress();
AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines->GetEntriesFast()));
- if(fLines->GetEntriesFast()>1){
- retcode = 0; // this last tracklet selection will be used
- delete fVert3D;
- fVert3D = 0;
- }
- else {
- retcode =1; // the previous tracklet selection will be used
- }
+ if(fLines->GetEntriesFast()>1) retcode=0; // this new tracklet selection is used
+ else retcode =1; // the previous tracklet selection will be used
}
else {
retcode = 0;
#ifndef ALIITSVERTEXER3D_H
#define ALIITSVERTEXER3D_H
-#include<AliITSVertexerZ.h>
+#include<AliITSVertexer.h>
///////////////////////////////////////////////////////////////////
// //
virtual ~AliITSVertexer3D();
virtual AliESDVertex* FindVertexForCurrentEvent(Int_t evnumb);
virtual void FindVertices();
- AliVertex *GetVertex3D() const {return fVert3D;}
+ AliVertex GetVertex3D() const {return fVert3D;}
virtual void MakeTracklet(Double_t *pA, Double_t *pB, Int_t &nolines);
virtual void MakeTracklet(Float_t *pA, Float_t *pB, Int_t &nolines);
virtual void PrintStatus() const;
AliITSVertexer3D(const AliITSVertexer3D& vtxr);
AliITSVertexer3D& operator=(const AliITSVertexer3D& /* vtxr */);
Int_t FindTracklets(Int_t evnumber, Int_t optCuts);
- void Find3DVertex();
Int_t Prepare3DVertex(Int_t optCuts);
+ void ResetVert3D();
TClonesArray *fLines; //! array of tracklets
- AliVertex *fVert3D; // 3D Vertex
+ AliVertex fVert3D; // 3D Vertex
Float_t fCoarseDiffPhiCut; // loose cut on DeltaPhi for RecPoint matching
Float_t fCoarseMaxRCut; // cut on tracklet DCA to Z axis
Float_t fMaxRCut; // cut on tracklet DCA to beam axis
}
totst /= nN;
Float_t cut = totst+totst2*2.;
+ // if(fDebug>1)cout<<"totst, totst2, cut: "<<totst<<", "<<totst2<<", "<<cut<<endl;
Float_t val1=hist->GetBinLowEdge(sepa);
Float_t val2=hist->GetBinLowEdge(1);
Float_t valm = 0.;
if((val2-val1)>deltaVal){
val1 = zmax-deltaVal/2.;
val2 = zmax+deltaVal/2.;
+ // if(fDebug>0)cout<<"val1 and val2 recomputed\n";
}
+ // if(fDebug>0)cout<<"Values for Z finding: "<<val1<<" "<<val2<<" "<<val2-val1<<endl;
fZFound=0;
fZsig=0;
nN=0;
if(firipixe>1){
rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
zave=zave/firipixe;
+ // if(fDebug>1)cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++\n Number of firing pixels: "<<firipixe<<endl;
}
else {
fZFound = -200;
Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.);
zlim2=zlim1 + sepa/10.;
TH1F *zvdis = new TH1F("z_ev","zv distr",sepa,zlim1,zlim2);
+ /*
+ if(fDebug>0){
+ cout<<"Z limits: "<<zlim1<<" "<<zlim2<<"; Bins= "<<sepa<<endl;
+ cout<<"Bin width: "<<zvdis->GetBinWidth(1)*10000.<<" micron\n";
+ }
+ */
Int_t sizarr=100;
TArrayF *zval = new TArrayF(sizarr);
// TArrayF *curv = new TArrayF(sizarr);
Int_t ncoinc=0;
for(Int_t module= fFirstL1; module<=fLastL1;module++){
+ // if(fDebug>0)cout<<"processing module "<<module<<" \r";
branch->GetEvent(module);
Int_t nrecp1 = itsRec->GetEntries();
TObjArray *poiL1 = new TObjArray(nrecp1);
Float_t r1=TMath::Sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
Float_t phi1 = TMath::ATan2(gc[1],gc[0]);
if(phi1<0)phi1=2*TMath::Pi()+phi1;
+ // if(fDebug>1)cout<<"module "<<module<<" "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<phi1<<" \n";
for(Int_t modul2=fFirstL2; modul2<=fLastL2; modul2++){
branch->GetEvent(modul2);
Int_t nrecp2 = itsRec->GetEntries();
}
delete poiL1;
} // loop on modules
+ /*
+ if(fDebug>0){
+ cout<<endl<<"Number of coincidences = "<<ncoinc<<endl;
+ }
+ */
// EvalZ(zvdis,sepa,ncoinc,zval,curv);
EvalZ(zvdis,sepa,ncoinc,zval);
delete zvdis;
WriteCurrentVertex();
}
else {
+ /*
+ if(fDebug>0){
+ cout<<"Vertex not found for event "<<i<<endl;
+ cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
+ }
+ */
}
}
}
#include "AliITSVertexerZ.h"
#include<TBranch.h>
#include<TClonesArray.h>
-#include<TFile.h>
#include<TH1.h>
#include <TString.h>
#include<TTree.h>
#include "AliITSDetTypeRec.h"
#include "AliITSRecPoint.h"
#include "AliITSZPoint.h"
-#include "AliHeader.h"
-#include "AliGenEventHeader.h"
/////////////////////////////////////////////////////////////////
// this class implements a fast method to determine
return fCurrentVertex;
}
-
-
-
//______________________________________________________________________
void AliITSVertexerZ::VertexZFinder(Int_t evnumber){
// Defines the AliESDVertex for the current event
TTree *tR = itsLoader->TreeR();
detTypeRec.SetTreeAddressR(tR);
TClonesArray *itsRec = 0;
- // lc and gc are local and global coordinates for layer 1
- Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
- Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
+ // lc1 and gc1 are local and global coordinates for layer 1
+ Float_t lc1[3]; for(Int_t ii=0; ii<3; ii++) lc1[ii]=0.;
+ Float_t gc1[3]; for(Int_t ii=0; ii<3; ii++) gc1[ii]=0.;
// lc2 and gc2 are local and global coordinates for layer 2
Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
Int_t nrpL1 = 0;
Int_t nrpL2 = 0;
+
// By default fFirstL1=0 and fLastL1=79
- // This loop counts the number of recpoints on layer1 (central modules)
for(Int_t module= fFirstL1; module<=fLastL1;module++){
- // Keep only central modules
- // if(module%4==0 || module%4==3)continue;
- // cout<<"Procesing module "<<module<<" ";
branch->GetEvent(module);
- // cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
nrpL1+= itsRec->GetEntries();
detTypeRec.ResetRecPoints();
}
//By default fFirstL2=80 and fLastL2=239
- //This loop counts the number of RP on layer 2
for(Int_t module= fFirstL2; module<=fLastL2;module++){
branch->GetEvent(module);
nrpL2+= itsRec->GetEntries();
fCurrentVertex = new AliESDVertex(0.,5.3,-2);
return;
}
- // The vertex finding is attempted only if the number of RP is !=0 on
- // both layers
- Float_t *xc1 = new Float_t [nrpL1]; // coordinates of the L1 Recpoints
- Float_t *yc1 = new Float_t [nrpL1];
- Float_t *zc1 = new Float_t [nrpL1];
- Float_t *phi1 = new Float_t [nrpL1];
- Float_t *err1 = new Float_t [nrpL1];
- Float_t *xc2 = new Float_t [nrpL2]; // coordinates of the L1 Recpoints
- Float_t *yc2 = new Float_t [nrpL2];
- Float_t *zc2 = new Float_t [nrpL2];
- Float_t *phi2 = new Float_t [nrpL2];
- Float_t *err2 = new Float_t [nrpL2];
- Int_t ind = 0;// running index for RP
// Force a coarse bin size of 200 microns if the number of clusters on layer 2
// is low
if(nrpL2<fPPsetting[0])SetBinWidthCoarse(fPPsetting[1]);
if(fZCombc)delete fZCombc;
fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
- // Loop on modules of layer 1
-
- for(Int_t module= fFirstL1; module<=fLastL1;module++){
- // if(module%4==0 || module%4==3)continue;
- branch->GetEvent(module);
- Int_t nrecp1 = itsRec->GetEntries();
- for(Int_t j=0;j<nrecp1;j++){
- AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
- // Local coordinates of this recpoint
- lc[0]=recp->GetDetLocalX();
- lc[2]=recp->GetDetLocalZ();
- geom->LtoG(module,lc,gc);
- // Global coordinates of this recpoints
- gc[0]-=fNominalPos[0]; // Possible beam offset in the bending plane
- gc[1]-=fNominalPos[1]; // " "
- xc1[ind]=gc[0];
- yc1[ind]=gc[1];
- zc1[ind]=gc[2];
- // azimuthal angle is computed in the interval 0 --> 2*pi
- phi1[ind] = TMath::ATan2(gc[1],gc[0]);
- if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
- err1[ind]=recp->GetSigmaZ2();
- ind++;
- }
- detTypeRec.ResetRecPoints();
- }
- ind = 0; // the running index is reset for Layer 2
- for(Int_t module= fFirstL2; module<=fLastL2;module++){
- branch->GetEvent(module);
- Int_t nrecp2 = itsRec->GetEntries();
- for(Int_t j=0;j<nrecp2;j++){
- AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
- lc[0]=recp->GetDetLocalX();
- lc[2]=recp->GetDetLocalZ();
- geom->LtoG(module,lc,gc);
- gc[0]-=fNominalPos[0];
- gc[1]-=fNominalPos[1];
- xc2[ind]=gc[0];
- yc2[ind]=gc[1];
- zc2[ind]=gc[2];
- phi2[ind] = TMath::ATan2(gc[1],gc[0]);
- if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
- err2[ind]=recp->GetSigmaZ2();
- ind++;
- }
- detTypeRec.ResetRecPoints();
- }
-
-/* Test the ffect of mutiple scatternig on error. Negligible
+ /* Test the ffect of mutiple scatternig on error. Negligible
// Multiple scattering
Float_t beta=1.,pmed=0.875; //pmed=875 MeV (for tracks with dphi<0.01 rad)
Float_t beta2=beta*beta;
Float_t theta2BP=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dBP);
Float_t theta2L1=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dL1);
*/
- TClonesArray *points = new TClonesArray("AliITSZPoint",nrpL1*nrpL2);
+ Int_t maxdim=TMath::Min(nrpL1*nrpL2,50000); // temporary; to limit the size in PbPb
+ TClonesArray *points = new TClonesArray("AliITSZPoint",maxdim);
TClonesArray &pts = *points;
Int_t nopoints =0;
- for(Int_t i=0;i<nrpL1;i++){ // loop on L1 RP
- Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]); // radius L1 RP
- for(Int_t j=0;j<nrpL2;j++){ // loop on L2 RP
- Float_t diff = TMath::Abs(phi2[j]-phi1[i]); // diff in azimuth
- if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; //diff<pi
- if(diff<fDiffPhiMax){ // cut on 10 milliradians by def.
- Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]); // radius L2 RP
-// Float_t tgth=(zc2[j]-zc1[i])/(r2-r1); // slope
- Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1); //Z @ null radius
- Float_t ezr0q=(r2*r2*err1[i]+r1*r1*err2[j])/(r2-r1)/(r2-r1); //error on Z @ null radius
- /*
+ for(Int_t modul1= fFirstL1; modul1<=fLastL1;modul1++){ // Loop on modules of layer 1
+ UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
+ branch->GetEvent(modul1);
+ Int_t nrecp1 = itsRec->GetEntries();
+ TClonesArray *prpl1 = new TClonesArray("AliITSRecPoint",nrecp1);
+ prpl1->SetOwner();
+ TClonesArray &rpl1 = *prpl1;
+ for(Int_t j=0;j<nrecp1;j++){
+ AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
+ new(rpl1[j])AliITSRecPoint(*recp);
+ }
+ detTypeRec.ResetRecPoints();
+ for(Int_t j1=0;j1<nrecp1;j1++){
+ AliITSRecPoint *recp = (AliITSRecPoint*)prpl1->At(j1);
+ lc1[0]=recp->GetDetLocalX();
+ lc1[2]=recp->GetDetLocalZ();
+ geom->LtoG(modul1,lc1,gc1);
+ // Global coordinates of this recpoints
+ gc1[0]-=fNominalPos[0]; // Possible beam offset in the bending plane
+ gc1[1]-=fNominalPos[1]; // " "
+ Float_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
+ Float_t phi1 = TMath::ATan2(gc1[1],gc1[0]);
+ if(phi1<0)phi1+=2*TMath::Pi();
+ Float_t zc1=gc1[2];
+ Float_t erz1=recp->GetSigmaZ2();
+ for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
+ for(Int_t k=0;k<4;k++){
+ Int_t ladmod=fLadders[ladder-1]+ladl2;
+ if(ladmod>geom->GetNladders(2)) ladmod=ladmod-geom->GetNladders(2);
+ Int_t modul2=geom->GetModuleIndex(2,ladmod,k+1);
+ branch->GetEvent(modul2);
+ Int_t nrecp2 = itsRec->GetEntries();
+ for(Int_t j2=0;j2<nrecp2;j2++){
+ recp = (AliITSRecPoint*)itsRec->At(j2);
+ lc2[0]=recp->GetDetLocalX();
+ lc2[2]=recp->GetDetLocalZ();
+ geom->LtoG(modul2,lc2,gc2);
+ gc2[0]-=fNominalPos[0];
+ gc2[1]-=fNominalPos[1];
+ Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
+ Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
+ if(phi2<0)phi2+=2*TMath::Pi();
+ Float_t zc2=gc2[2];
+ Float_t erz2=recp->GetSigmaZ2();
+
+ Float_t diff = TMath::Abs(phi2-phi1);
+ if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
+ if(diff<fDiffPhiMax){
+ // Float_t tgth=(zc2[j]-zc1[i])/(r2-r1); // slope (used for multiple scattering)
+ Float_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius
+ Float_t ezr0q=(r2*r2*erz1+r1*r1*erz2)/(r2-r1)/(r2-r1); //error on Z @ null radius
+ /*
+ // Multiple scattering
ezr0q+=r1*r1*(1+tgth*tgth)*theta2L1/2; // multiple scattering in layer 1
ezr0q+=rBP*rBP*(1+tgth*tgth)*theta2BP/2; // multiple scattering in beam pipe
*/
- new(pts[nopoints++])AliITSZPoint(zr0,ezr0q);
-
- fZCombc->Fill(zr0);
+ if(nopoints<maxdim) new(pts[nopoints++])AliITSZPoint(zr0,ezr0q);
+ fZCombc->Fill(zr0);
+ }
+ }
+ detTypeRec.ResetRecPoints();
+ }
}
}
+ delete prpl1;
}
- delete [] xc1;
- delete [] yc1;
- delete [] zc1;
- delete [] phi1;
- delete [] err1;
- delete [] xc2;
- delete [] yc2;
- delete [] zc2;
- delete [] phi2;
- delete [] err2;
points->Sort();
} while(niter<10 && TMath::Abs((zm-lim1)-(lim2-zm))>fTolerance);
fCurrentVertex = new AliESDVertex(zm,ezm,ncontr);
fCurrentVertex->SetTitle("vertexer: B");
+ delete points;
ResetHistograms();
itsLoader->UnloadRecPoints();
return;
}
+
+
//_____________________________________________________________________
void AliITSVertexerZ::ResetHistograms(){
// delete TH1 data members
virtual void FindVertices();
virtual void PrintStatus() const;
void SetDiffPhiMax(Float_t pm = 0.01){fDiffPhiMax = pm;}
- void SetVtxStart(Float_t x, Float_t y){
- fX0=x; fY0=y;
- }
void ConfigIterations(Int_t noiter=3,Float_t *ptr=0);
void SetFirstLayerModules(Int_t m1 = 0, Int_t m2 = 79){fFirstL1 = m1; fLastL1 = m2;}
void SetSecondLayerModules(Int_t m1 = 80, Int_t m2 = 239){fFirstL2 = m1; fLastL2 = m2;}
Int_t fFirstL2; // first module of the second pixel layer used
Int_t fLastL2; // last module of the second pixel layer used
Float_t fDiffPhiMax; // Maximum delta phi allowed among corr. pixels
- Float_t fX0; // Nominal x coordinate of the vertex
- Float_t fY0; // Nominal y coordinate of the vertex
Float_t fZFound; //! found value for the current event
Float_t fZsig; //! RMS of Z
TH1F *fZCombc; //! histogram with coarse z distribution
Float_t fPhiDiffIter[5]; // Delta phi used in iterations
Float_t fWindowWidth; // Z window width for symmetrization
- ClassDef(AliITSVertexerZ,6);
+ ClassDef(AliITSVertexerZ,7);
};
#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff