X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSMeanVertexer.cxx;h=caa5e5927b48ebe68c9cb8129f98b944bebe3a2f;hb=ff7f662866547dcccdf38b3b952d324436791f61;hp=69548b9000dd88e845584850261fdba75e5507ec;hpb=308c2f7c71668bd4c7193564d1da2ebd974ea86c;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSMeanVertexer.cxx b/ITS/AliITSMeanVertexer.cxx index 69548b9000d..caa5e5927b4 100644 --- a/ITS/AliITSMeanVertexer.cxx +++ b/ITS/AliITSMeanVertexer.cxx @@ -1,9 +1,12 @@ #include +#include +#include +#include #include "AliGeomManager.h" -#include "AliHeader.h" #include "AliITSDetTypeRec.h" #include "AliITSInitGeometry.h" #include "AliITSMeanVertexer.h" +#include "AliITSRecPointContainer.h" #include "AliITSLoader.h" #include "AliLog.h" #include "AliRawReader.h" @@ -11,15 +14,13 @@ #include "AliRawReaderRoot.h" #include "AliRunLoader.h" #include "AliITSVertexer3D.h" +#include "AliITSVertexer3DTapan.h" #include "AliESDVertex.h" #include "AliMeanVertex.h" -#include "AliMultiplicity.h" - -const TString AliITSMeanVertexer::fgkMVFileNameDefault = "MeanVertex.root"; - +//#include "AliCodeTimer.h" +const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000; ClassImp(AliITSMeanVertexer) - /////////////////////////////////////////////////////////////////////// // // // Class to compute vertex position using SPD local reconstruction // @@ -27,277 +28,313 @@ ClassImp(AliITSMeanVertexer) // 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():TObject(), -fLoaderFileName(), -fGeometryFileName(), -fMVFileName(), -fRawReader(), -fRunLoader(), +AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(), +fDetTypeRec(NULL), +fVertexXY(NULL), +fVertexZ(NULL), fNoEventsContr(0), -fTotTracklets(0.), -fAverTracklets(0.), -fSigmaOnAverTracks(0.), +fTotContributors(0.), +fAverContributors(0.), fFilterOnContributors(0), -fFilterOnTracklets(0), -fWriteVertices(kFALSE) +fMode(mode), +fVertexer(NULL), +fAccEvents(fgkMaxNumOfEvents), +fVertArray("AliESDVertex",fgkMaxNumOfEvents), +fClu0(NULL), +fIndex(0), +fErrXCut(0.), +fRCut(0.), +fZCutDiamond(0.), +fLowSPD0(0), +fHighSPD0(0), +fMultH(NULL), +fErrXH(NULL), +fMultHa(NULL), +fErrXHa(NULL), +fDistH(NULL), +fContrH(NULL), +fContrHa(NULL) { // Default Constructor - for(Int_t i=0;i<3;i++){ - fWeighPos[i] = 0.; - fWeighSig[i] = 0.; - fAverPos[i] = 0.; - for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.; - } - + SetZFiducialRegion(); // sets fZCutDiamond to its default value + Reset(kFALSE,kFALSE); + + // Histograms initialization + 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); + fVertexXY->SetXTitle("X , cm"); + fVertexXY->SetYTitle("Y , cm"); + fVertexXY->SetOption("colz"); + 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;iGetXaxis()->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(); + } } - SetLoaderFileName(); - SetGeometryFileName(); - - Init(filename); -} -//______________________________________________________________________ -AliITSMeanVertexer::AliITSMeanVertexer(TString &filename, - TString &loaderfilename, - TString &geometryfilename):TObject(), -fLoaderFileName(), -fGeometryFileName(), -fMVFileName(fgkMVFileNameDefault), -fRawReader(), -fRunLoader(), -fNoEventsContr(0), -fTotTracklets(0.), -fAverTracklets(0.), -fSigmaOnAverTracks(0.), -fFilterOnContributors(0), -fFilterOnTracklets(0), -fWriteVertices(kTRUE) -{ - // Standard constructor with explicit geometry file name assignment for(Int_t i=0;i<3;i++){ + fWeighPosSum[i] = 0.; + fWeighSigSum[i] = 0.; + fAverPosSum[i] = 0.; fWeighPos[i] = 0.; fWeighSig[i] = 0.; fAverPos[i] = 0.; for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.; + for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.; + fNoEventsContr=0; + fTotContributors = 0.; } - SetLoaderFileName(loaderfilename); - SetGeometryFileName(geometryfilename); - Init(filename); } - //______________________________________________________________________ -void AliITSMeanVertexer::Init(TString &filename){ - // Initialization part common to different constructors - if(filename.IsNull()){ - AliFatal("Please, provide a valid file name for raw data file\n"); - } - // if file name ends with root a raw reader ROOT is assumed - if(filename.EndsWith(".root")){ - fRawReader = new AliRawReaderRoot(filename); - } - else { // DATE raw reader is assumed - fRawReader = new AliRawReaderDate(filename); - fRawReader->SelectEvents(7); - } - fRunLoader = AliRunLoader::Open(fLoaderFileName.Data(),AliConfig::GetDefaultEventFolderName(),"recreate"); - fRunLoader->MakeTree("E"); - Int_t iEvent = 0; - while (fRawReader->NextEvent()) { - fRunLoader->SetEventNumber(iEvent); - fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(), - iEvent, iEvent); - fRunLoader->MakeTree("H"); - fRunLoader->TreeE()->Fill(); - iEvent++; - } - fRawReader->RewindEvents(); - fRunLoader->SetNumberOfEventsPerFile(iEvent); - fRunLoader->WriteHeader("OVERWRITE"); - Int_t retval = AliConfig::Instance()->AddDetector(fRunLoader->GetEventFolder(),"ITS","ITS"); - if(retval != 0)AliFatal("Not able to add ITS detector"); - AliITSLoader *loader = new AliITSLoader("ITS",fRunLoader->GetEventFolder()->GetName()); - fRunLoader->AddLoader(loader); - fRunLoader->CdGAFile(); - fRunLoader->Write(0, TObject::kOverwrite); +Bool_t AliITSMeanVertexer::Init() { + // Initialize filters // Initialize geometry + // Initialize ITS classes - AliGeomManager::LoadGeometry(fGeometryFileName.Data()); + AliGeomManager::LoadGeometry(); + if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE; AliITSInitGeometry initgeom; AliITSgeom *geom = initgeom.CreateAliITSgeom(); + if (!geom) return kFALSE; printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data()); - loader->SetITSgeom(geom); - // Initialize filter values to their defaults - SetFilterOnContributors(); - SetFilterOnTracklets(); -} -//______________________________________________________________________ -AliITSMeanVertexer::AliITSMeanVertexer(const AliITSMeanVertexer &vtxr) : TObject(vtxr), -fLoaderFileName(vtxr.fLoaderFileName), -fGeometryFileName(vtxr.fGeometryFileName), -fMVFileName(vtxr.fMVFileName), -fRawReader(vtxr.fRawReader), -fRunLoader(vtxr.fRunLoader), -fNoEventsContr(vtxr.fNoEventsContr), -fTotTracklets(vtxr.fTotTracklets), -fAverTracklets(vtxr.fAverTracklets), -fSigmaOnAverTracks(vtxr.fSigmaOnAverTracks), -fFilterOnContributors(vtxr.fFilterOnContributors), -fFilterOnTracklets(vtxr.fFilterOnTracklets), -fWriteVertices(vtxr.fWriteVertices) -{ - // Copy constructor - // Copies are not allowed. The method is protected to avoid misuse. - AliFatal("Copy constructor not allowed\n"); + fDetTypeRec = new AliITSDetTypeRec(); + fDetTypeRec->SetLoadOnlySPDCalib(kTRUE); + fDetTypeRec->SetITSgeom(geom); + fDetTypeRec->SetDefaults(); + fDetTypeRec->SetDefaultClusterFindersV2(kTRUE); -} + // Initialize filter values to their defaults + SetFilterOnContributors(); + SetCutOnErrX(); + SetCutOnR(); + SetCutOnCls(); -//______________________________________________________________________ -AliITSMeanVertexer& AliITSMeanVertexer::operator=(const AliITSMeanVertexer& /* vtxr */ ){ - // Assignment operator - AliError("Assignment operator not allowed\n"); - return *this; + // Instatiate vertexer + if (!fMode) { + fVertexer = new AliITSVertexer3DTapan(1000); + } + else { + fVertexer = new AliITSVertexer3D(); + fVertexer->SetDetTypeRec(fDetTypeRec); + AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer; + alias->SetWideFiducialRegion(fZCutDiamond,0.5); + alias->SetNarrowFiducialRegion(0.5,0.5); + alias->SetDeltaPhiCuts(0.5,0.025); + alias->SetDCACut(0.1); + alias->SetPileupAlgo(3); + fVertexer->SetComputeMultiplicity(kFALSE); + } + return kTRUE; } //______________________________________________________________________ AliITSMeanVertexer::~AliITSMeanVertexer() { // Destructor - delete fRawReader; - delete fRunLoader; - + delete fDetTypeRec; + delete fVertexXY; + delete fVertexZ; + delete fMultH; + delete fErrXH; + delete fMultHa; + delete fErrXHa; + delete fDistH; + delete fContrH; + delete fContrHa; + delete fVertexer; } //______________________________________________________________________ -void AliITSMeanVertexer::Reconstruct(){ +Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){ // Performs SPD local reconstruction - AliITSLoader* loader = static_cast(fRunLoader->GetLoader("ITSLoader")); - if (!loader) { - AliFatal("ITS loader not found"); - return; - } - AliITSDetTypeRec* rec = new AliITSDetTypeRec(); - rec->SetITSgeom(loader->GetITSgeom()); - rec->SetDefaults(); + // and vertex finding + // returns true in case a vertex is found - rec->SetDefaultClusterFindersV2(kTRUE); + // 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"); - Int_t nEvents = fRunLoader->GetNumberOfEvents(); - fRawReader->RewindEvents(); - for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) { - fRawReader->NextEvent(); - fRunLoader->GetEvent(iEvent); - AliDebug(1,Form(">>>>>>> Processing event number: %d",iEvent)); - loader->LoadRecPoints("update"); - loader->CleanRecPoints(); - loader->MakeRecPointsContainer(); - TTree *tR = loader->TreeR(); - if(!tR){ - AliFatal("Tree R pointer not found - Abort \n"); - break; + Bool_t vtxOK = kFALSE; + UInt_t nl1=0; + AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree); + if (!fMode) { + if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE; + } + else { + 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(); } - rec->DigitsToRecPoints(fRawReader,tR,"SPD"); - rec->ResetRecPoints(); - rec->ResetClusters(); - loader->WriteRecPoints("OVERWRITE"); - loader->UnloadRecPoints(); + 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; + return vtxOK; } //______________________________________________________________________ -void AliITSMeanVertexer::DoVertices(){ - // Loop on all events and compute 3D vertices - AliITSLoader* loader = static_cast(fRunLoader->GetLoader("ITSLoader")); - AliITSVertexer3D *dovert = new AliITSVertexer3D(); - AliESDVertex *vert = 0; - Int_t nevents = fRunLoader->TreeE()->GetEntries(); - for(Int_t i=0; iGetEvent(i); - TTree* cltree = loader->TreeR(); - vert = dovert->FindVertexForCurrentEvent(cltree); - AliMultiplicity *mult = dovert->GetMultiplicity(); - if(Filter(vert,mult)){ - AddToMean(vert); - if(fWriteVertices){ - loader->PostVertex(vert); - loader->WriteVertices(); - } - } - else { - if(vert)delete vert; +void AliITSMeanVertexer::WriteVertices(const char *filename){ + // Compute mean vertex and + // store it along with the histograms + // in a file + + TFile fmv(filename,"update"); + Bool_t itisOK = kFALSE; + if(ComputeMean(kFALSE)){ + if(ComputeMean(kTRUE)){ + Double_t cov[6]; + cov[0] = fAverPosSq[0][0]; // variance x + cov[1] = fAverPosSq[0][1]; // cov xy + cov[2] = fAverPosSq[1][1]; // variance y + cov[3] = fAverPosSq[0][2]; // cov xz + cov[4] = fAverPosSq[1][2]; // cov yz + cov[5] = fAverPosSq[2][2]; // variance z + // 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(fAverPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos"); + + mv.Write(mv.GetName(),TObject::kOverwrite); + vtx.Write(vtx.GetName(),TObject::kOverwrite); + itisOK = kTRUE; } } - if(ComputeMean()){ - Double_t cov[6]; - cov[0] = fAverPosSq[0][0]; // variance x - cov[1] = fAverPosSq[0][1]; // cov xy - cov[2] = fAverPosSq[1][1]; // variance y - 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,nevents,fTotTracklets,fAverTracklets,fSigmaOnAverTracks); - mv.SetTitle("Mean Vertex"); - mv.SetName("Meanvertex"); - AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets())); - AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks())); - TFile fmv(fMVFileName.Data(),"recreate"); - mv.Write(); - fmv.Close(); - } - else { + + if(!itisOK){ AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr)); } - delete dovert; + + 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(); } //______________________________________________________________________ -Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){ +Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){ // Apply selection criteria to events Bool_t status = kFALSE; - if(!vert || !mult)return status; + if(!vert)return status; + // Remove vertices reconstructed with vertexerZ + if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status; Int_t ncontr = vert->GetNContributors(); - Int_t ntracklets = mult->GetNumberOfTracklets(); AliDebug(1,Form("Number of contributors = %d",ncontr)); - AliDebug(1,Form("Number of tracklets = %d",ntracklets)); - if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE; - fAverTracklets += ntracklets; - fSigmaOnAverTracks += ntracklets*ntracklets; + Double_t ex = vert->GetXRes(); + fMultH->Fill(mult); + fErrXH->Fill(ex*1000.); + fContrH->Fill((Float_t)ncontr); + if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && multFill(mult); + fErrXHa->Fill(ex*1000.); + fContrHa->Fill((Float_t)ncontr); + } return status; } @@ -311,39 +348,86 @@ void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){ for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE; if(!goon)return; for(Int_t i=0;i<3;i++){ - fWeighPos[i]+=currentPos[i]/currentSigma[i]/currentSigma[i]; - fWeighSig[i]+=1./currentSigma[i]/currentSigma[i]; - fAverPos[i]+=currentPos[i]; + fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i]; + fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i]; + fAverPosSum[i]+=currentPos[i]; } for(Int_t i=0;i<3;i++){ for(Int_t j=i;j<3;j++){ - fAverPosSq[i][j] += currentPos[i] * currentPos[j]; + fAverPosSqSum[i][j] += currentPos[i] * currentPos[j]; } } + + fTotContributors+=vert->GetNContributors(); + fVertexXY->Fill(currentPos[0],currentPos[1]); + fVertexZ->Fill(currentPos[2]); + fNoEventsContr++; } //______________________________________________________________________ -Bool_t AliITSMeanVertexer::ComputeMean(){ - // compute mean vertex - if(fNoEventsContr < 2) return kFALSE; - Double_t nevents = fNoEventsContr; - for(Int_t i=0;i<3;i++){ - fWeighPos[i] /= fWeighSig[i]; - fWeighSig[i] = 1./TMath::Sqrt(fWeighSig[i]); - fAverPos[i] /= nevents; - } - for(Int_t i=0;i<3;i++){ - for(Int_t j=i;j<3;j++){ - fAverPosSq[i][j] /= (nevents -1.); - fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j]; - } - } - fTotTracklets = fAverTracklets; // total number of tracklets used - fAverTracklets /= nevents; - fSigmaOnAverTracks /= (nevents - 1); - Double_t tmp = nevents/(nevents -1.)*fAverTracklets*fAverTracklets; - fSigmaOnAverTracks -= tmp; - fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks); - return kTRUE; -} + Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){ + // 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]; + + Int_t istart = oldnumbevt; + if(killOutliers)istart = 0; + if(killOutliers && preliminary){ + preliminary = kFALSE; + Reset(kTRUE,kFALSE); + } + oldnumbevt = fVertArray.GetEntries(); + for(Int_t i=istart;iGetXv()-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); + } + + if(!fAccEvents.TestBitNumber(i))continue; + vert=(AliESDVertex*)fVertArray[i]; + AddToMean(vert); + } + 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]; + fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]); + fAverPos[i] = fAverPosSum[i]/nevents; + } + for(Int_t i=0;i<3;i++){ + for(Int_t j=i;j<3;j++){ + fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.); + fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j]; + } + } + if(killOutliers)ResetArray(); + fAverContributors = fTotContributors/nevents; + return kTRUE; + }