#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
+#include <TTree.h>
#include "AliGeomManager.h"
#include "AliITSDetTypeRec.h"
#include "AliITSInitGeometry.h"
#include "AliITSVertexer3DTapan.h"
#include "AliESDVertex.h"
#include "AliMeanVertex.h"
+//#include "AliCodeTimer.h"
const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
ClassImp(AliITSMeanVertexer)
-
///////////////////////////////////////////////////////////////////////
// //
// Class to compute vertex position using SPD local reconstruction //
// is built and saved //
// Individual vertices can be optionally stored //
// Origin: M.Masera (masera@to.infn.it) //
-// Usage:
-// AliITSMeanVertexer mv("RawDataFileName");
-// mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
-// .... optional setters ....
-// mv.Reconstruct(); // SPD local reconstruction
-// mv.DoVertices();
-// Resulting AliMeanVertex object is stored in a file named fMVFileName
+// Usage: //
+// Class used by the ITSSPDVertexDiamondda.cxx detector algorithm //
///////////////////////////////////////////////////////////////////////
-
/* $Id$ */
-
//______________________________________________________________________
AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
fDetTypeRec(NULL),
fIndex(0),
fErrXCut(0.),
fRCut(0.),
+fZCutDiamond(0.),
fLowSPD0(0),
-fHighSPD0(0)
+fHighSPD0(0),
+fMultH(NULL),
+fErrXH(NULL),
+fMultHa(NULL),
+fErrXHa(NULL),
+fDistH(NULL),
+fContrH(NULL),
+fContrHa(NULL)
{
// Default Constructor
- Reset();
+ SetZFiducialRegion(); // sets fZCutDiamond to its default value
+ Reset(kFALSE,kFALSE);
// Histograms initialization
- const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
- const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
+ const Float_t xLimit = 0.6, yLimit = 0.6, zLimit = 50.0;
+ const Float_t xDelta = 0.003, yDelta = 0.003, zDelta = 0.2;
fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
fVertexZ->SetXTitle("Z , cm");
+ fErrXH = new TH1F("errorX","Error X - before cuts",100,0.,99.);
+ fMultH = new TH1F("multh","mult on layer 1 before cuts",2400,0.,7200.);
+ fErrXHa = new TH1F("errorXa","Error X - after Filter",100,0.,99.);
+ fErrXHa->SetLineColor(kRed);
+ fMultHa = new TH1F("multha","mult on layer 1 - after Filter",2400,0.,7200.);
+ fMultHa->SetLineColor(kRed);
+ fDistH = new TH1F("disth","distance (mm)",100,0.,30.);
+ fContrH = new TH1F("contrib","Number of contributors - before cuts",300,0.5,300.5);
+ fContrHa = new TH1F("contriba","Number of contributors - after cuts",300,0.5,300.5);
+ fContrHa->SetLineColor(kRed);
+
fClu0 = new UInt_t [fgkMaxNumOfEvents];
for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
fAccEvents.ResetAllBits(kTRUE);
}
//______________________________________________________________________
-void AliITSMeanVertexer::Reset(){
+void AliITSMeanVertexer::Reset(Bool_t redefine2D,Bool_t complete){
// Reset averages
+ // Both arguments are kFALSE when the method is called by the constructor
+ // redefine2D is TRUE when the method is called by ComputeMean at the second
+ // pass.
+ // redefine2D is FALSE and complete is TRUE when the method is called
+ // after a failure cof ComputeMean(kTRUE)
+
+ static Int_t counter = 0;
+ if(fVertexXY && redefine2D){
+ counter++;
+ Double_t rangeX=fVertexXY->GetXaxis()->GetXmax()-fVertexXY->GetXaxis()->GetXmin();
+ Double_t rangeY=fVertexXY->GetYaxis()->GetXmax()-fVertexXY->GetYaxis()->GetXmin();
+ Int_t nx=fVertexXY->GetNbinsX();
+ Int_t ny=fVertexXY->GetNbinsY();
+ delete fVertexXY;
+ Double_t xmi=fWeighPos[0]-rangeX/2.;
+ Double_t xma=fWeighPos[0]+rangeX/2.;
+ Double_t ymi=fWeighPos[1]-rangeY/2.;
+ Double_t yma=fWeighPos[1]+rangeY/2.;
+ fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",nx,xmi,xma,ny,ymi,yma);
+ fVertexXY->SetXTitle("X , cm");
+ fVertexXY->SetYTitle("Y , cm");
+ fVertexXY->SetOption("colz");
+ fVertexXY->SetDirectory(0);
+ }
+ else if(fVertexXY && !redefine2D){
+ fVertexXY->Reset();
+ }
+ if(fVertexZ){
+ fVertexZ->Reset();
+ fDistH->Reset();
+ if(complete){
+ fErrXH->Reset();
+ fMultH->Reset();
+ fErrXHa->Reset();
+ fMultHa->Reset();
+ fContrH->Reset();
+ fContrHa->Reset();
+ }
+ }
for(Int_t i=0;i<3;i++){
fWeighPosSum[i] = 0.;
fWeighSigSum[i] = 0.;
for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
fNoEventsContr=0;
fTotContributors = 0.;
- if(fVertexXY)fVertexXY->Reset();
- if(fVertexZ)fVertexZ->Reset();
}
}
//______________________________________________________________________
fVertexer = new AliITSVertexer3D();
fVertexer->SetDetTypeRec(fDetTypeRec);
AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
- alias->SetWideFiducialRegion(40.,0.5);
+ alias->SetWideFiducialRegion(fZCutDiamond,0.5);
alias->SetNarrowFiducialRegion(0.5,0.5);
alias->SetDeltaPhiCuts(0.5,0.025);
alias->SetDCACut(0.1);
delete fDetTypeRec;
delete fVertexXY;
delete fVertexZ;
+ delete fMultH;
+ delete fErrXH;
+ delete fMultHa;
+ delete fErrXHa;
+ delete fDistH;
+ delete fContrH;
+ delete fContrHa;
delete fVertexer;
}
// returns true in case a vertex is found
// Run SPD cluster finder
+ static Int_t evcount = -1;
+ if(evcount <0){
+ evcount++;
+ AliInfo(Form("Low and high cuts on SPD L0 clusters %d , %d \n",fLowSPD0,fHighSPD0));
+ AliInfo(Form("Reconstruct: cut on errX %f \n",fErrXCut));
+ }
+// AliCodeTimerAuto("",0);
AliITSRecPointContainer::Instance()->PrepareToRead();
TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
nl1=rpcont->GetNClustersInLayerFast(1);
if(Filter(vtx,nl1)) vtxOK = kTRUE;
+ /* if(vtx){
+ if(vtxOK){
+ printf("The vertex is OK\n");
+ }
+ else {
+ printf("The vertex is NOT OK\n");
+ }
+ vtx->PrintStatus();
+ }
+ else {
+ printf("The vertex was not reconstructed\n");
+ } */
}
delete clustersTree;
if (vtxOK && fMode){
new(fVertArray[fIndex])AliESDVertex(*vtx);
fClu0[fIndex]=nl1;
+ fAccEvents.SetBitNumber(fIndex);
fIndex++;
+ if(fIndex>=fgkMaxNumOfEvents){
+ if(ComputeMean(kFALSE)){
+ if(ComputeMean(kTRUE))AliInfo("Mean vertex computed");
+ }
+ }
}
if (vtx) delete vtx;
cov[3] = fAverPosSq[0][2]; // cov xz
cov[4] = fAverPosSq[1][2]; // cov yz
cov[5] = fAverPosSq[2][2]; // variance z
- AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
+ // We use standard average and not weighed averag now
+ // AliMeanVertex is apparently not taken by the preprocessor; only
+ // the AliESDVertex object is retrieved
+ // AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
+ AliMeanVertex mv(fAverPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
mv.SetTitle("Mean Vertex");
mv.SetName("MeanVertex");
// we have to add chi2 here
- AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
+ // AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
+ AliESDVertex vtx(fAverPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
mv.Write(mv.GetName(),TObject::kOverwrite);
vtx.Write(vtx.GetName(),TObject::kOverwrite);
fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
+ fMultH->Write(fMultH->GetName(),TObject::kOverwrite);
+ fErrXH->Write(fErrXH->GetName(),TObject::kOverwrite);
+ fMultHa->Write(fMultHa->GetName(),TObject::kOverwrite);
+ fErrXHa->Write(fErrXHa->GetName(),TObject::kOverwrite);
+ fDistH->Write(fDistH->GetName(),TObject::kOverwrite);
+ fContrH->Write(fContrH->GetName(),TObject::kOverwrite);
+ fContrHa->Write(fContrHa->GetName(),TObject::kOverwrite);
fmv.Close();
}
Int_t ncontr = vert->GetNContributors();
AliDebug(1,Form("Number of contributors = %d",ncontr));
Double_t ex = vert->GetXRes();
+ fMultH->Fill(mult);
+ fErrXH->Fill(ex*1000.);
+ fContrH->Fill((Float_t)ncontr);
if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) &&
((ex*1000.)<fErrXCut)) {
status = kTRUE;
- fAccEvents.SetBitNumber(fIndex);
+ fMultHa->Fill(mult);
+ fErrXHa->Fill(ex*1000.);
+ fContrHa->Fill((Float_t)ncontr);
}
return status;
}
// compute mean vertex
// called once with killOutliers = kFALSE and then again with
// killOutliers = kTRUE
+
+ static Bool_t preliminary = kTRUE;
+ static Int_t oldnumbevt = 0;
+ if(!(preliminary || killOutliers))return kTRUE; //ComputeMean is
+ // executed only once with argument kFALSE
Double_t wpos[3];
for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
- Reset();
- for(Int_t i=0;i<fVertArray.GetEntries();i++){
+
+ Int_t istart = oldnumbevt;
+ if(killOutliers)istart = 0;
+ if(killOutliers && preliminary){
+ preliminary = kFALSE;
+ Reset(kTRUE,kFALSE);
+ }
+ oldnumbevt = fVertArray.GetEntries();
+ for(Int_t i=istart;i<fVertArray.GetEntries();i++){
AliESDVertex* vert = NULL;
// second pass
Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
dist=sqrt(dist)*10.; // distance in mm
+ fDistH->Fill(dist);
if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
}
vert=(AliESDVertex*)fVertArray[i];
AddToMean(vert);
}
- if(fNoEventsContr < 2) return kFALSE;
+ Bool_t bad = ((!killOutliers) && fNoEventsContr < 5) || (killOutliers && fNoEventsContr <2);
+ if(bad) {
+ if(killOutliers){
+// when the control reachs this point, the preliminary evaluation of the
+// diamond region has to be redone. So a general reset must be issued
+// if this occurs, it is likely that the cuts are badly defined
+ ResetArray();
+ Reset(kFALSE,kTRUE);
+ preliminary = kTRUE;
+ oldnumbevt = 0;
+ }
+ return kFALSE;
+ }
Double_t nevents = fNoEventsContr;
for(Int_t i=0;i<3;i++){
fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
}
}
-
+ if(killOutliers)ResetArray();
fAverContributors = fTotContributors/nevents;
return kTRUE;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff