X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSMeanVertexer.cxx;h=388549e26e06fb05957027d0f51d331efa8287e7;hb=6d64cce3e519d7e242b33e116dbbb4b2ef539127;hp=69548b9000dd88e845584850261fdba75e5507ec;hpb=308c2f7c71668bd4c7193564d1da2ebd974ea86c;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSMeanVertexer.cxx b/ITS/AliITSMeanVertexer.cxx index 69548b9000d..388549e26e0 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,13 +14,11 @@ #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"; - +const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000; ClassImp(AliITSMeanVertexer) /////////////////////////////////////////////////////////////////////// @@ -39,265 +40,202 @@ ClassImp(AliITSMeanVertexer) /* $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.), +fLowSPD0(0), +fHighSPD0(0) { // 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.; - } - -} - -//______________________________________________________________________ -AliITSMeanVertexer::AliITSMeanVertexer(TString &filename):TObject(), -fLoaderFileName(), -fGeometryFileName(), -fMVFileName(fgkMVFileNameDefault), -fRawReader(), -fRunLoader(), -fNoEventsContr(0), -fTotTracklets(0.), -fAverTracklets(0.), -fSigmaOnAverTracks(0.), -fFilterOnContributors(0), -fFilterOnTracklets(0), -fWriteVertices(kTRUE) -{ - // Standard constructor + Reset(); - 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.; - } - SetLoaderFileName(); - SetGeometryFileName(); - - Init(filename); + // 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; + 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"); + fClu0 = new UInt_t [fgkMaxNumOfEvents]; + for(Int_t i=0;iReset(); + if(fVertexZ)fVertexZ->Reset(); } - 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(40.,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 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 + 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; - } - rec->DigitsToRecPoints(fRawReader,tR,"SPD"); - rec->ResetRecPoints(); - rec->ResetClusters(); - loader->WriteRecPoints("OVERWRITE"); - loader->UnloadRecPoints(); + 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; + } + delete clustersTree; + if (vtxOK && fMode){ + new(fVertArray[fIndex])AliESDVertex(*vtx); + fClu0[fIndex]=nl1; + fIndex++; + } + 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 + AliMeanVertex mv(fWeighPos,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"); + + 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); + 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(); + if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && multGetNContributors(); + 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 + Double_t wpos[3]; + for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i]; + Reset(); + for(Int_t i=0;iGetXv()-wpos[0])*(vert->GetXv()-wpos[0]); + dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]); + dist=sqrt(dist)*10.; // distance in mm + if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE); + } + + if(!fAccEvents.TestBitNumber(i))continue; + vert=(AliESDVertex*)fVertArray[i]; + AddToMean(vert); + } + if(fNoEventsContr < 2) 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]; + } + } + + fAverContributors = fTotContributors/nevents; + return kTRUE; + }