From e0ea2e34437e22acf620e342388fd64548496bb1 Mon Sep 17 00:00:00 2001 From: dainese Date: Mon, 25 Jul 2011 16:05:27 +0000 Subject: [PATCH] Class and macro to update the OCDB with the MeanVertex during pass0 (Davide) --- PWG1/CalibMacros/Pass0/MeanVertexCalibPass0.C | 19 + PWG1/ITS/AliMeanVertexPreprocessorOffline.cxx | 419 ++++++++++++++++++ PWG1/ITS/AliMeanVertexPreprocessorOffline.h | 32 ++ PWG1/PWG1LinkDef.h | 1 + 4 files changed, 471 insertions(+) create mode 100644 PWG1/CalibMacros/Pass0/MeanVertexCalibPass0.C create mode 100644 PWG1/ITS/AliMeanVertexPreprocessorOffline.cxx create mode 100644 PWG1/ITS/AliMeanVertexPreprocessorOffline.h diff --git a/PWG1/CalibMacros/Pass0/MeanVertexCalibPass0.C b/PWG1/CalibMacros/Pass0/MeanVertexCalibPass0.C new file mode 100644 index 00000000000..4e3ff6eb75a --- /dev/null +++ b/PWG1/CalibMacros/Pass0/MeanVertexCalibPass0.C @@ -0,0 +1,19 @@ +//___________________________________________________________________________ + +LoadLibraries() +{ + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + +} + +//___________________________________________________________________________ + +MakeOCDB(const Char_t *filename = "AliESDfriends_v1.root", const Char_t *dbString = "raw://alice/data/2011/OCDB/", Int_t runNB) +{ + LoadLibraries(); + AliMeanVertexPreprocessorOffline meanVertexCalib; + meanVertexCalib.ProcessOutput(filename, dbString, runNB); +} + + diff --git a/PWG1/ITS/AliMeanVertexPreprocessorOffline.cxx b/PWG1/ITS/AliMeanVertexPreprocessorOffline.cxx new file mode 100644 index 00000000000..f266e39bd70 --- /dev/null +++ b/PWG1/ITS/AliMeanVertexPreprocessorOffline.cxx @@ -0,0 +1,419 @@ +/************************************************************************** + * Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* + +*/ +// Mean Vertex preprocessor: +// 2) takes data after pass0 , +// processes it, and stores either to OCDB . +// +// Davide Caffarri + +#include "AliMeanVertexPreprocessorOffline.h" + +#include "AliCDBStorage.h" +#include "AliCDBMetaData.h" +#include "AliCDBManager.h" +#include "AliCDBEntry.h" +#include "AliLog.h" + +#include +#include +#include +#include +#include "TClass.h" + +#include "AliESDVertex.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TF1.h" +#include "TProfile.h" + + +ClassImp(AliMeanVertexPreprocessorOffline) + +//____________________________________________________ +AliMeanVertexPreprocessorOffline::AliMeanVertexPreprocessorOffline(): +TNamed("AliMeanVertexPreprocessorOffline","AliMeanVertexPreprocessorOffline") + +{ + //constructor +} +//____________________________________________________ + +AliMeanVertexPreprocessorOffline::~AliMeanVertexPreprocessorOffline() +{ + //destructor + +} +//____________________________________________________ +void AliMeanVertexPreprocessorOffline::ProcessOutput(const char *filename, const char *dbString, Int_t runNb){ + + TFile *file = TFile::Open(filename); + if (!file || !file->IsOpen()){ + AliError(Form("cannot open outputfile %s", filename)); + return; + } + + TList *list = (TList*)file->Get("MeanVertex"); + if (!list) { + AliError(Form("cannot find the list of histograms")); + return; + } + + // check on the histos + + TH1F *histTRKvtxX = (TH1F*)list->FindObject("hTRKVertexX"); + TH1F *histTRKvtxY = (TH1F*)list->FindObject("hTRKVertexY"); + TH1F *histTRKvtxZ = (TH1F*)list->FindObject("hTRKVertexZ"); + + TH1F *histSPDvtxX = (TH1F*)list->FindObject("hSPDVertexX"); + TH1F *histSPDvtxY = (TH1F*)list->FindObject("hSPDVertexY"); + TH1F *histSPDvtxZ = (TH1F*)list->FindObject("hSPDVertexZ"); + + Bool_t useTRKvtx = kTRUE; + Bool_t useITSSAvtx = kFALSE; + Bool_t useSPDvtx = kFALSE; + + if (!histTRKvtxX || !histTRKvtxY || !histTRKvtxZ){ + + useTRKvtx = kFALSE; + useITSSAvtx = kTRUE; + + histTRKvtxX = (TH1F*)list->FindObject("hITSSAVertexX"); + histTRKvtxY = (TH1F*)list->FindObject("hITSSAVertexY"); + histTRKvtxZ = (TH1F*)list->FindObject("hITSSAVertexZ"); + + + if (!histTRKvtxX || !histTRKvtxY || !histTRKvtxZ){ + + useITSSAvtx=kFALSE; + useSPDvtx=kTRUE; + + if (!histSPDvtxX || !histSPDvtxY || !histSPDvtxZ){ + AliError(Form("cannot find any histograms available")); + return; + } + } + } + + + if (useTRKvtx){ + + Float_t nEntriesX = histTRKvtxX->GetEffectiveEntries(); + Float_t nEntriesY = histTRKvtxY->GetEffectiveEntries(); + Float_t nEntriesZ = histTRKvtxZ->GetEffectiveEntries(); + + if (nEntriesX < 50. || nEntriesY<50. || nEntriesZ<50.) { + AliWarning(Form("TRK vertex histograms have too few entries for fitting")); + + if (useITSSAvtx){ + if (nEntriesX < 50. || nEntriesY<50. || nEntriesZ<50.) { + AliWarning(Form("ITSSA vertex histograms have too few entries for fitting")); + + if (useSPDvtx){ + + nEntriesX = histSPDvtxX->GetEffectiveEntries(); + nEntriesY = histSPDvtxY->GetEffectiveEntries(); + nEntriesZ = histSPDvtxZ->GetEffectiveEntries(); + + if (nEntriesX < 50. || nEntriesY<50. || nEntriesZ<50.) { + AliError(Form("Also SPD vertex histograms have too few entries for fitting, return")); + return; + } + } + } + } + } + } + + + Double_t xMeanVtx=0., yMeanVtx=0., zMeanVtx=0.; + Double_t xSigmaVtx=0., ySigmaVtx=0., zSigmaVtx=0.; + + TF1 *fitVtxX, *fitVtxY, *fitVtxZ; + + if (useTRKvtx || useITSSAvtx){ + histTRKvtxX ->Fit("gaus", "M"); + fitVtxX = histTRKvtxX -> GetFunction("gaus"); + xMeanVtx = fitVtxX -> GetParameter(1); + if (TMath::Abs(xMeanVtx) > 2.) { + xMeanVtx = 0.; + } + + histTRKvtxY ->Fit("gaus", "M"); + fitVtxY = histTRKvtxY -> GetFunction("gaus"); + yMeanVtx = fitVtxY -> GetParameter(1); + if (TMath::Abs(yMeanVtx) > 2.) { + yMeanVtx = 0.; + } + + histTRKvtxZ ->Fit("gaus", "M"); + fitVtxZ = histTRKvtxZ -> GetFunction("gaus"); + zMeanVtx = fitVtxZ -> GetParameter(1); + zSigmaVtx = fitVtxZ -> GetParameter(2); + if (TMath::Abs(zMeanVtx) > 8.) { + zMeanVtx = 0.; + } + + } + + + //check fits: compare histo mean with fit mean value + Double_t xHistoMean, yHistoMean, zHistoMean; + Double_t xHistoRMS, yHistoRMS, zHistoRMS; + + if (useTRKvtx || useITSSAvtx){ + xHistoMean = histTRKvtxX -> GetMean(); + xHistoRMS = histTRKvtxX ->GetRMS(); + + if ((TMath::Abs(xHistoMean-xMeanVtx) > 0.5)){ + AliWarning(Form("Possible problems with the fit mean very different from histo mean... using SPD vertex")); + useTRKvtx = kFALSE; + useITSSAvtx = kFALSE; + useSPDvtx = kTRUE; + } + + yHistoMean = histTRKvtxY -> GetMean(); + yHistoRMS = histTRKvtxY ->GetRMS(); + + if ((TMath::Abs(yHistoMean-yMeanVtx) > 0.5)){ + AliWarning(Form("Possible problems with the fit mean very different from histo mean... using SPD vertex")); + useTRKvtx = kFALSE; + useITSSAvtx = kFALSE; + useSPDvtx = kTRUE; + } + + zHistoMean = histTRKvtxZ -> GetMean(); + zHistoRMS = histTRKvtxZ ->GetRMS(); + + if ((TMath::Abs(zHistoMean-zMeanVtx) > 1.)){ + AliWarning(Form("Possible problems with the fit mean very different from histo mean... using SPD vertex")); + useTRKvtx = kFALSE; + useITSSAvtx = kFALSE; + useSPDvtx = kTRUE; + } + } + + + if (useSPDvtx){ + + histSPDvtxX ->Fit("gaus", "M"); + fitVtxX = histSPDvtxX -> GetFunction("gaus"); + xMeanVtx = fitVtxX -> GetParameter(1); + if (TMath::Abs(xMeanVtx) > 2.) { + xMeanVtx = 0.; + } + + histSPDvtxY ->Fit("gaus", "M"); + fitVtxY = histSPDvtxY -> GetFunction("gaus"); + yMeanVtx = fitVtxY -> GetParameter(1); + if (TMath::Abs(yMeanVtx) > 2.) { + yMeanVtx = 0.; + } + + histSPDvtxZ ->Fit("gaus", "M"); + fitVtxZ = histSPDvtxZ -> GetFunction("gaus"); + zMeanVtx = fitVtxZ -> GetParameter(1); + zSigmaVtx = fitVtxZ -> GetParameter(2); + if (TMath::Abs(zMeanVtx) > 7.) { + zMeanVtx = 0.; + zSigmaVtx = 5.; + } + + } + + + + //check with online position + + if (useTRKvtx || useITSSAvtx){ + AliCDBManager *manCheck = AliCDBManager::Instance(); + manCheck->SetDefaultStorage("raw://"); + manCheck->SetRun(runNb); + + AliCDBEntry *entr = manCheck->Get("GRP/Calib/MeanVertexSPD"); + AliESDVertex *vtxOnline = (AliESDVertex*)entr->GetObject(); + + Double_t posOnline[3], sigmaOnline[3]; + posOnline[0] = vtxOnline->GetX(); + posOnline[1] = vtxOnline->GetY(); + posOnline[2] = vtxOnline->GetZ(); + + sigmaOnline[0] = vtxOnline->GetXRes(); + sigmaOnline[1] = vtxOnline->GetYRes(); + sigmaOnline[2] = vtxOnline->GetZRes(); + + //vtxOnline->GetSigmaXYZ(sigmaOnline); + + if ((TMath::Abs(posOnline[0]-xMeanVtx) > 0.1) || (TMath::Abs(posOnline[1]-yMeanVtx) > 0.1) || (TMath::Abs(posOnline[2]-zMeanVtx) > 1.)){ + AliWarning(Form("vertex offline far from the online one")); + return; + } + + } + + Float_t meanMult = 40.; + Float_t p2 = 1.4; + Float_t resolVtx = 0.05; + + Double_t xSigmaMult, ySigmaMult, corrXZ, corrYZ, lumiRegSquaredX, lumiRegSquaredY; + Double_t covarXZ=0., covarYZ=0.; + Bool_t histosLumiReg = kTRUE; + + TF1 *sigmaFitX, *sigmaFitY, *corrFit; + + TH1F *histTRKdefMultX=0; + TH1F *histTRKdefMultY=0; + TH2F *histTRKVertexXZ=0; + TH2F *histTRKVertexYZ=0; + + if ((!histTRKdefMultX) || (!histTRKdefMultY)){ + AliWarning(Form("histos for lumi reg calculation not found, default value setted")); + xSigmaMult=0.0120; + ySigmaMult=0.0120; + histosLumiReg=kFALSE; + + } + + if (histosLumiReg){ + if (useTRKvtx){ + + histTRKdefMultX = (TH1F*)list->FindObject("hTRKVertexXdefMult"); + histTRKdefMultY = (TH1F*)list->FindObject("hTRKVertexYdefMult"); + histTRKVertexXZ = (TH2F*)list->FindObject("hTRKVertexXZ"); + histTRKVertexYZ = (TH2F*)list->FindObject("hTRKVertexYZ"); + + } + + if (useITSSAvtx){ + + histTRKdefMultX = (TH1F*)list->FindObject("hITSSAVertexXdefMult"); + histTRKdefMultY = (TH1F*)list->FindObject("hITSSAVertexYdefMult"); + histTRKVertexXZ = (TH2F*)list->FindObject("hITSSAVertexXZ"); + histTRKVertexYZ = (TH2F*)list->FindObject("hITSSAVertexYZ"); + } + + + + histTRKdefMultX -> Fit("gaus", "M", "", -1., 1.); + sigmaFitX = histTRKdefMultX -> GetFunction("gaus"); + xSigmaMult = sigmaFitX->GetParameter(2); + + + lumiRegSquaredX = (xSigmaMult*xSigmaMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2))); + if (lumiRegSquaredX < 0) { + AliWarning(Form("Problems with luminosiy region determination, update of the postion only")); + xSigmaMult = 0.; + xSigmaVtx = 0.0120; + } + + if (lumiRegSquaredX > 0 && lumiRegSquaredX < 0.0005){ + xSigmaVtx = TMath::Sqrt(lumiRegSquaredX); + xSigmaVtx = xSigmaVtx*1.2; + } + + + histTRKdefMultY -> Fit("gaus", "M"); + sigmaFitY = histTRKdefMultY -> GetFunction("gaus"); + ySigmaMult = sigmaFitY->GetParameter(2); + + lumiRegSquaredY= (ySigmaMult*ySigmaMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2))); + if (lumiRegSquaredY < 0) { + AliWarning(Form("Problems with luminosiy region determination, update of the postion only")); + ySigmaMult = 0.; + ySigmaVtx = 0.0120; + } + + if (lumiRegSquaredY > 0 && lumiRegSquaredY < 0.0005){ + ySigmaVtx = TMath::Sqrt(lumiRegSquaredY); + ySigmaVtx = ySigmaVtx*1.2; + } + + TProfile *htrkXZ = histTRKVertexXZ ->ProfileY(); + htrkXZ -> Fit("pol1", "M", "", -10., 10.); + corrFit = htrkXZ->GetFunction("pol1"); + corrXZ = corrFit->GetParameter(1); + + if (TMath::Abs(corrXZ) > 0.01) { + AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix")); + corrXZ =0.; + } + else{ + covarXZ = corrXZ * zSigmaVtx*zSigmaVtx; + + } + + TProfile *htrkYZ = histTRKVertexYZ ->ProfileY(); + htrkYZ -> Fit("pol1", "M", "", -10., 10.); + corrFit = htrkYZ->GetFunction("pol1"); + corrYZ = corrFit->GetParameter(1); + + if (TMath::Abs(corrYZ) > 0.01) { + AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix")); + corrYZ =0.; + } + else{ + covarYZ = corrYZ*zSigmaVtx*zSigmaVtx; + } + + } + + Double_t position[3], covMatrix[6]; + Double_t chi2=1.; + Int_t nContr=1; + + position[0] = xMeanVtx; + position[1] = yMeanVtx; + position[2] = zMeanVtx; + + covMatrix[0] = xSigmaVtx; + covMatrix[1] = 0.; //xy + covMatrix[2] = ySigmaVtx; + covMatrix[3] = covarXZ; + covMatrix[4] = covarYZ; + covMatrix[5] = zSigmaVtx; + + + Printf ("%f, %f, %f, %f", xSigmaVtx, ySigmaVtx, covarXZ, covarYZ); + + AliESDVertex *vertex = new AliESDVertex(position, covMatrix, chi2, nContr, "vertex"); + + AliCDBManager *cdb = AliCDBManager::Instance(); + AliCDBStorage *sto = cdb->GetStorage(dbString); + if (!sto) { + AliError(Form("cannot get storage %s", dbString)); + return; + } + + AliCDBId id("GRP/Calib/MeanVertex", runNb, runNb); + + AliCDBMetaData metaData; + metaData.SetBeamPeriod(0); //check!!!! + metaData.SetResponsible("Davide Caffarri"); + metaData.SetComment("Mean Vertex object used in reconstruction"); + + if (!cdb->Put(vertex, id, &metaData)) { + AliError(Form("error while putting object in storage %s", dbString)); + } + + delete vertex; + + + +} + + diff --git a/PWG1/ITS/AliMeanVertexPreprocessorOffline.h b/PWG1/ITS/AliMeanVertexPreprocessorOffline.h new file mode 100644 index 00000000000..cae1a314d7d --- /dev/null +++ b/PWG1/ITS/AliMeanVertexPreprocessorOffline.h @@ -0,0 +1,32 @@ +#ifndef ALI_MEANVERTEX_PREPROCESSOROFFLINE_H +#define ALI_MEANVERTEX_PREPRECESSOROFFLINE_H + +/* Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id: AliMeanVertexPreprocessor.h $ */ + + +// Mean vertex preprocessor. +// Davide Caffarri +// +#include "TNamed.h" + +class AliMeanVertexPreprocessorOffline: public TNamed +{ + public: + AliMeanVertexPreprocessorOffline(); + virtual ~AliMeanVertexPreprocessorOffline(); + + void ProcessOutput(const char *filename, const char *dbString, Int_t runNb); + + + private: + AliMeanVertexPreprocessorOffline(const AliMeanVertexPreprocessorOffline & proc); // copy constructor + AliMeanVertexPreprocessorOffline& operator=(const AliMeanVertexPreprocessorOffline&); //operator + + + ClassDef(AliMeanVertexPreprocessorOffline, 1); +}; + +#endif diff --git a/PWG1/PWG1LinkDef.h b/PWG1/PWG1LinkDef.h index 52180e5a29e..a971e22d9d6 100644 --- a/PWG1/PWG1LinkDef.h +++ b/PWG1/PWG1LinkDef.h @@ -72,6 +72,7 @@ #pragma link C++ class AliSPDUtils+; #pragma link C++ class AliAnalysisTaskdEdxSSDQA+; #pragma link C++ class AliMeanVertexCalibTask+; +#pragma link C++ class AliMeanVertexPreprocessorOffline+; #pragma link C++ class AliRelAlignerKalmanArray+; #pragma link C++ class AliAnalysisTaskITSTPCalignment+; -- 2.43.0