[u/mrichter/AliRoot.git] / PWGPP / ITS / AliMeanVertexPreprocessorOffline.cxx

/**************************************************************************
 * 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 <TTimeStamp.h>
#include <TFile.h>
#include <TObjString.h>
#include <TNamed.h>
#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; 
	}

    if (!dbString){
		AliError(Form("no OCDB path found, return"));
	    return;
	}
    
	
	TList *list = (TList*)file->Get("MeanVertex");				
	
	TH1F *histTRKvtxX = 0x0;
	TH1F *histTRKvtxY = 0x0;
	TH1F *histTRKvtxZ = 0x0;
	
	TH1F *histSPDvtxX = 0x0;
	TH1F *histSPDvtxY = 0x0;
	TH1F *histSPDvtxZ = 0x0;
	
	
	Bool_t useTRKvtx = kTRUE;
	Bool_t useITSSAvtx = kFALSE;
	Bool_t useSPDvtx = kFALSE;
	Bool_t spdAvailable = kTRUE;
	Bool_t writeMeanVertexSPD = kFALSE;
	Bool_t vertexerSPD3Doff=kFALSE;
	
	
    if (!list) {
		
		histTRKvtxX = (TH1F*)file->Get("hTRKVertexX");
		histTRKvtxY = (TH1F*)file->Get("hTRKVertexY");
		histTRKvtxZ = (TH1F*)file->Get("hTRKVertexZ");
		
		histSPDvtxX = (TH1F*)file->Get("hSPDVertexX");
		histSPDvtxY = (TH1F*)file->Get("hSPDVertexY");
		histSPDvtxZ = (TH1F*)file->Get("hSPDVertexZ");
		
		if (!histTRKvtxX || !histTRKvtxY || !histTRKvtxZ){
			
			useTRKvtx = kFALSE;
			useITSSAvtx = kTRUE;
			
			histTRKvtxX = (TH1F*)file->FindObject("hITSSAVertexX");
			histTRKvtxY = (TH1F*)file->FindObject("hITSSAVertexY");
			histTRKvtxZ = (TH1F*)file->FindObject("hITSSAVertexZ");
		
			if (!histTRKvtxX || !histTRKvtxY || !histTRKvtxZ){
				
				useITSSAvtx=kFALSE;
				useSPDvtx=kTRUE;
				
				if (!histSPDvtxX || !histSPDvtxY || !histSPDvtxZ){
					AliError(Form("cannot find any histograms available from file"));
					return;
				}	
			}
		}
	}	
	
	else{
		
		histTRKvtxX = (TH1F*)list->FindObject("hTRKVertexX");
		histTRKvtxY = (TH1F*)list->FindObject("hTRKVertexY");
		histTRKvtxZ = (TH1F*)list->FindObject("hTRKVertexZ");
		
		histSPDvtxX = (TH1F*)list->FindObject("hSPDVertexX");
		histSPDvtxY = (TH1F*)list->FindObject("hSPDVertexY");
		histSPDvtxZ = (TH1F*)list->FindObject("hSPDVertexZ");
	
		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 from list"));
					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.) {
			AliError(Form("TRK vertex histograms have too few entries for fitting"));
			useTRKvtx=kFALSE;
			useSPDvtx = kTRUE;	
		}
	}
	if (useITSSAvtx){
	  
		Float_t nEntriesX = histTRKvtxX->GetEffectiveEntries(); 					 
		Float_t nEntriesY = histTRKvtxY->GetEffectiveEntries(); 			 
		Float_t nEntriesZ = histTRKvtxZ->GetEffectiveEntries(); 			 
	
		if (nEntriesX < 50. || nEntriesY<50. || nEntriesZ<50.) {
			AliError(Form("ITSSA vertex histograms have too few entries for fitting"));
			useITSSAvtx=kFALSE;
			useSPDvtx=kTRUE;
		}
	}
		
	Float_t nEntriesX = histSPDvtxX->GetEffectiveEntries(); 					 
	Float_t nEntriesY = histSPDvtxY->GetEffectiveEntries(); 			 
	Float_t nEntriesZ = histSPDvtxZ->GetEffectiveEntries(); 

	if (nEntriesX < 50. || nEntriesY<50. || nEntriesZ<50.) {
		spdAvailable = kFALSE;
		if ((useTRKvtx==kFALSE) && (useITSSAvtx==kFALSE)){
			AliError(Form("Also SPD vertex histograms have too few entries for fitting, return"));
			return;		 
		}
	}
	
	if((nEntriesX == 0.)&&(nEntriesY==0.) && (nEntriesZ>0.)){
		vertexerSPD3Doff = kTRUE;
		AliWarning("Vertexer SPD 3D off");
	}
	
	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.;
			writeMeanVertexSPD=kTRUE;
		}	

		histTRKvtxY ->Fit("gaus", "M");
		fitVtxY = histTRKvtxY -> GetFunction("gaus");
		yMeanVtx = fitVtxY -> GetParameter(1);
		if (TMath::Abs(yMeanVtx) > 2.) {
			yMeanVtx = 0.;
			writeMeanVertexSPD=kTRUE;
		}	
		
		histTRKvtxZ ->Fit("gaus", "M");
		fitVtxZ = histTRKvtxZ -> GetFunction("gaus");
		zMeanVtx = fitVtxZ -> GetParameter(1);
		zSigmaVtx = fitVtxZ -> GetParameter(2);
		if ((TMath::Abs(zMeanVtx) > 20.) || (zSigmaVtx>12.)) {
		  zMeanVtx = histTRKvtxZ->GetMean();
		  zSigmaVtx = histTRKvtxZ->GetRMS();
		  writeMeanVertexSPD=kTRUE;
		}	
		
	}
	
	
	//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) && (spdAvailable) && (!vertexerSPD3Doff)){
	  
	  histSPDvtxX ->Fit("gaus", "M");
	  fitVtxX = histSPDvtxX -> GetFunction("gaus");
	  xMeanVtx = fitVtxX -> GetParameter(1);
	  if (TMath::Abs(xMeanVtx) > 2.) {
			xMeanVtx = 0.;
			writeMeanVertexSPD=kTRUE;
	  }
		
	  histSPDvtxY ->Fit("gaus", "M");
	  fitVtxY = histSPDvtxY -> GetFunction("gaus");
	  yMeanVtx = fitVtxY -> GetParameter(1);
		if (TMath::Abs(yMeanVtx) > 2.) {
			yMeanVtx = 0.;
			writeMeanVertexSPD=kTRUE;
		}	
		
		histSPDvtxZ ->Fit("gaus", "M");
		fitVtxZ = histSPDvtxZ -> GetFunction("gaus");
		zMeanVtx = fitVtxZ -> GetParameter(1);
		zSigmaVtx = fitVtxZ -> GetParameter(2);
		if ((TMath::Abs(zMeanVtx) > 20.) || (zSigmaVtx>12.)) {
		  zMeanVtx = histSPDvtxZ ->GetMean();
		  zSigmaVtx = histSPDvtxZ->GetRMS();
		  writeMeanVertexSPD = kTRUE;
		}	
				
	}
	else if ((useSPDvtx) && (!spdAvailable)){
		AliError(Form("Difference between trkVtx and online one, SPD histos not enough entry or SPD 3D vertex off. Writing Mean Vertex SPD"));
		writeMeanVertexSPD = kTRUE;	
	}
	
	
	//check with online position
	
	 Double_t posOnline[3], sigmaOnline[3];
	
	if (useTRKvtx || useITSSAvtx || writeMeanVertexSPD){
		AliCDBManager *manCheck = AliCDBManager::Instance();
		manCheck->SetDefaultStorage("raw://");
		manCheck->SetRun(runNb);
		
		AliCDBEntry *entr = manCheck->Get("GRP/Calib/MeanVertexSPD");
		if(entr) {
		  AliESDVertex *vtxOnline = (AliESDVertex*)entr->GetObject();
		
		  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"));
		  }
		}
	}
	
	
	if (writeMeanVertexSPD){
		
		AliWarning(Form("Writing Mean Vertex SPD, Mean Vertex not available"));
				   
		Double_t sigma[3]={0.0150, 0.0150, zSigmaVtx};
			 
		AliESDVertex  *vertex =  new AliESDVertex(posOnline, sigma, "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 (!sto->Put(vertex, id, &metaData)) {
			AliError(Form("Error while putting object in storage %s", dbString));
		}
		
		delete vertex;
		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 highMultEnvironment = kTRUE;
	
	TF1 *sigmaFitX, *sigmaFitY, *corrFit;
		
	TH1F *histTRKdefMultX=0;
	TH1F *histTRKdefMultY=0;
	TH1F *histTRKHighMultX=0;
	TH1F *histTRKHighMultY=0;
	TH2F *histTRKVertexXZ=0;
	TH2F *histTRKVertexYZ=0;
	
	
	if (useTRKvtx){
	  if (list){
		  histTRKdefMultX = (TH1F*)list->FindObject("hTRKVertexXdefMult");
		  histTRKdefMultY = (TH1F*)list->FindObject("hTRKVertexYdefMult");
		  histTRKHighMultX = (TH1F*)list->FindObject("hTRKVertexXHighMult");
		  histTRKHighMultY = (TH1F*)list->FindObject("hTRKVertexYHighMult");
		  histTRKVertexXZ = (TH2F*)list->FindObject("hTRKVertexXZ");
		  histTRKVertexYZ = (TH2F*)list->FindObject("hTRKVertexYZ");
	  }
		
	  else {
		  histTRKdefMultX = (TH1F*)file->Get("hTRKVertexXdefMult");
		  histTRKdefMultY = (TH1F*)file->Get("hTRKVertexYdefMult");
		  histTRKHighMultX = (TH1F*)file->Get("hTRKVertexXHighMult");
		  histTRKHighMultY = (TH1F*)file->Get("hTRKVertexYHighMult");
		  histTRKVertexXZ = (TH2F*)file->Get("hTRKVertexXZ");
		  histTRKVertexYZ = (TH2F*)file->Get("hTRKVertexYZ");
	  }
		
	}	
	
	if (useITSSAvtx){
	  if (list){
		  histTRKdefMultX = (TH1F*)list->FindObject("hITSSAVertexXdefMult");
		  histTRKdefMultY = (TH1F*)list->FindObject("hITSSAVertexYdefMult");
		  histTRKHighMultX = (TH1F*)list->FindObject("hITSSAVertexXHighMult");
		  histTRKHighMultY = (TH1F*)list->FindObject("hITSSAVertexYHighMult");
		  histTRKVertexXZ = (TH2F*)list->FindObject("hITSSAVertexXZ");
		  histTRKVertexYZ = (TH2F*)list->FindObject("hITSSAVertexYZ");
	  }
	  
	  else {
		  histTRKdefMultX = (TH1F*)file->Get("hITSSAVertexXdefMult");
		  histTRKdefMultY = (TH1F*)file->Get("hITSSAVertexYdefMult");
		  histTRKHighMultX = (TH1F*)file->Get("hITSSAVertexXHighMult");
		  histTRKHighMultY = (TH1F*)file->Get("hITSSAVertexYHighMult");
		  histTRKVertexXZ = (TH2F*)file->Get("hITSSAVertexXZ");
		  histTRKVertexYZ = (TH2F*)file->Get("hITSSAVertexYZ");
	  }
	}
		
	
	if ((histTRKHighMultX) && (histTRKHighMultY)){
		
		Float_t nEntriesHighMultX = histTRKHighMultX->GetEffectiveEntries();
		Float_t nEntriesHighMultY = histTRKHighMultY->GetEffectiveEntries();
		
		if ((nEntriesHighMultX >100) && (nEntriesHighMultY>100)) {
			AliWarning(Form("Setting High Mulitplicity environment"));	
	     	highMultEnvironment = kTRUE;
		}
		else {
			AliWarning(Form("Setting Low Mulitplicity environment"));
			highMultEnvironment=kFALSE;			
		}
	}
	
	else{
		AliWarning(Form("No histos found, setting Low Mulitplicity environment"));
		highMultEnvironment=kFALSE;
		
	}
	
	if (highMultEnvironment==kFALSE){
		
		if ((!histTRKdefMultX) || (!histTRKdefMultY) || (histTRKdefMultX->GetEntries() < 40.) || (histTRKdefMultY->GetEntries() < 40.)){
			AliWarning(Form("histos for lumi reg calculation not found, default value setted"));
			xSigmaVtx=0.0120;
			ySigmaVtx=0.0120;
		} else {
	
			histTRKdefMultX -> Fit("gaus", "M", "", -0.4, 0.4);
			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", "", -0.2, 0.6);
			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;
			}
			
		}
	}
		
	if (highMultEnvironment==kTRUE){
	
		histTRKHighMultX -> Fit("gaus", "M", "", -0.4, 0.4);
		sigmaFitX = histTRKHighMultX -> GetFunction("gaus");
		xSigmaMult = sigmaFitX->GetParameter(2);
		
		if ((xSigmaMult <0) || (xSigmaMult>0.03)){
			AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
			xSigmaMult = 0.;
			xSigmaVtx = 0.0120;
		}
        else{
			xSigmaVtx = xSigmaMult;
		    xSigmaVtx = xSigmaVtx*1.2;
		}
			
		histTRKHighMultY -> Fit("gaus", "M", "", -0.2, 0.5);
		sigmaFitY = histTRKHighMultY -> GetFunction("gaus");
		ySigmaMult = sigmaFitY->GetParameter(2);
		
		if ((ySigmaMult <0) || (ySigmaMult>0.03)){
			AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
			ySigmaMult = 0.;
			ySigmaVtx = 0.0120;
		}
        else{
			ySigmaVtx = ySigmaMult;
		    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*xSigmaVtx;
  covMatrix[1] = 0.; //xy
  covMatrix[2] = ySigmaVtx*ySigmaVtx;
  covMatrix[3] = covarXZ;
  covMatrix[4] = covarYZ;
  covMatrix[5] = zSigmaVtx*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 (!sto->Put(vertex, id, &metaData)) {
	AliError(Form("Error while putting object in storage %s", dbString));
   }
 
  delete vertex;
 
  
}
Commit	Line	Data
e0ea2e34	1	/**************************************************************************
	2	* Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/*
	17
	18	*/
	19	// Mean Vertex preprocessor:
	20	// 2) takes data after pass0 ,
	21	// processes it, and stores either to OCDB .
	22	//
	23	// Davide Caffarri
	24
	25	#include "AliMeanVertexPreprocessorOffline.h"
	26
	27	#include "AliCDBStorage.h"
	28	#include "AliCDBMetaData.h"
	29	#include "AliCDBManager.h"
	30	#include "AliCDBEntry.h"
	31	#include "AliLog.h"
	32
	33	#include <TTimeStamp.h>
	34	#include <TFile.h>
	35	#include <TObjString.h>
	36	#include <TNamed.h>
	37	#include "TClass.h"
	38
	39	#include "AliESDVertex.h"
	40	#include "TH1F.h"
	41	#include "TH2F.h"
	42	#include "TF1.h"
	43	#include "TProfile.h"
	44
	45
	46	ClassImp(AliMeanVertexPreprocessorOffline)
	47
	48	//____________________________________________________
	49	AliMeanVertexPreprocessorOffline::AliMeanVertexPreprocessorOffline():
	50	TNamed("AliMeanVertexPreprocessorOffline","AliMeanVertexPreprocessorOffline")
	51
	52	{
	53	//constructor
	54	}
	55	//____________________________________________________
	56
	57	AliMeanVertexPreprocessorOffline::~AliMeanVertexPreprocessorOffline()
	58	{
	59	//destructor
	60
	61	}
	62	//____________________________________________________
	63	void AliMeanVertexPreprocessorOffline::ProcessOutput(const char filename, const char dbString, Int_t runNb){
	64
65	TFile *file = TFile::Open(filename);
66	if (!file \|\| !file->IsOpen()){
67	AliError(Form("cannot open outputfile %s", filename));
68	return;
69	}
70
792b5897	71	if (!dbString){
	72	AliError(Form("no OCDB path found, return"));
	73	return;
e0ea2e34	74	}
792b5897	75
	76
	77	TList list = (TList)file->Get("MeanVertex");
	78
	79	TH1F *histTRKvtxX = 0x0;
	80	TH1F *histTRKvtxY = 0x0;
	81	TH1F *histTRKvtxZ = 0x0;
	82
	83	TH1F *histSPDvtxX = 0x0;
	84	TH1F *histSPDvtxY = 0x0;
	85	TH1F *histSPDvtxZ = 0x0;
	86
e0ea2e34	87
	88	Bool_t useTRKvtx = kTRUE;
	89	Bool_t useITSSAvtx = kFALSE;
	90	Bool_t useSPDvtx = kFALSE;
be5ba8df	91	Bool_t spdAvailable = kTRUE;
be5ba8df	92	Bool_t writeMeanVertexSPD = kFALSE;
a680415a	93	Bool_t vertexerSPD3Doff=kFALSE;
e0ea2e34	94
792b5897	95
792b5897	96	if (!list) {
e0ea2e34	97
792b5897	98	histTRKvtxX = (TH1F*)file->Get("hTRKVertexX");
	99	histTRKvtxY = (TH1F*)file->Get("hTRKVertexY");
	100	histTRKvtxZ = (TH1F*)file->Get("hTRKVertexZ");
e0ea2e34	101
792b5897	102	histSPDvtxX = (TH1F*)file->Get("hSPDVertexX");
	103	histSPDvtxY = (TH1F*)file->Get("hSPDVertexY");
	104	histSPDvtxZ = (TH1F*)file->Get("hSPDVertexZ");
	105
	106	if (!histTRKvtxX \|\| !histTRKvtxY \|\| !histTRKvtxZ){
	107
	108	useTRKvtx = kFALSE;
	109	useITSSAvtx = kTRUE;
	110
	111	histTRKvtxX = (TH1F*)file->FindObject("hITSSAVertexX");
	112	histTRKvtxY = (TH1F*)file->FindObject("hITSSAVertexY");
	113	histTRKvtxZ = (TH1F*)file->FindObject("hITSSAVertexZ");
e0ea2e34	114
792b5897	115	if (!histTRKvtxX \|\| !histTRKvtxY \|\| !histTRKvtxZ){
	116
	117	useITSSAvtx=kFALSE;
	118	useSPDvtx=kTRUE;
	119
	120	if (!histSPDvtxX \|\| !histSPDvtxY \|\| !histSPDvtxZ){
	121	AliError(Form("cannot find any histograms available from file"));
	122	return;
	123	}
	124	}
	125	}
	126	}
	127
	128	else{
e0ea2e34	129
792b5897	130	histTRKvtxX = (TH1F*)list->FindObject("hTRKVertexX");
	131	histTRKvtxY = (TH1F*)list->FindObject("hTRKVertexY");
	132	histTRKvtxZ = (TH1F*)list->FindObject("hTRKVertexZ");
e0ea2e34	133
792b5897	134	histSPDvtxX = (TH1F*)list->FindObject("hSPDVertexX");
	135	histSPDvtxY = (TH1F*)list->FindObject("hSPDVertexY");
	136	histSPDvtxZ = (TH1F*)list->FindObject("hSPDVertexZ");
	137
	138	if (!histTRKvtxX \|\| !histTRKvtxY \|\| !histTRKvtxZ){
	139
	140	useTRKvtx = kFALSE;
	141	useITSSAvtx = kTRUE;
e0ea2e34	142
792b5897	143	histTRKvtxX = (TH1F*)list->FindObject("hITSSAVertexX");
	144	histTRKvtxY = (TH1F*)list->FindObject("hITSSAVertexY");
	145	histTRKvtxZ = (TH1F*)list->FindObject("hITSSAVertexZ");
	146
	147	if (!histTRKvtxX \|\| !histTRKvtxY \|\| !histTRKvtxZ){
	148
	149	useITSSAvtx=kFALSE;
	150	useSPDvtx=kTRUE;
	151
	152	if (!histSPDvtxX \|\| !histSPDvtxY \|\| !histSPDvtxZ){
	153	AliError(Form("cannot find any histograms available from list"));
	154	return;
	155	}
	156	}
	157
	158	}
	159	}
e0ea2e34	160
	161
	162	if (useTRKvtx){
adca66d2	163
792b5897	164	Float_t nEntriesX = histTRKvtxX->GetEffectiveEntries();
e0ea2e34	165	Float_t nEntriesY = histTRKvtxY->GetEffectiveEntries();
	166	Float_t nEntriesZ = histTRKvtxZ->GetEffectiveEntries();
	167
	168	if (nEntriesX < 50. \|\| nEntriesY<50. \|\| nEntriesZ<50.) {
adca66d2	169	AliError(Form("TRK vertex histograms have too few entries for fitting"));
ed15b00c	170	useTRKvtx=kFALSE;
ed15b00c	171	useSPDvtx = kTRUE;
adca66d2	172	}
	173	}
	174	if (useITSSAvtx){
	175
792b5897	176	Float_t nEntriesX = histTRKvtxX->GetEffectiveEntries();
adca66d2	177	Float_t nEntriesY = histTRKvtxY->GetEffectiveEntries();
	178	Float_t nEntriesZ = histTRKvtxZ->GetEffectiveEntries();
	179
	180	if (nEntriesX < 50. \|\| nEntriesY<50. \|\| nEntriesZ<50.) {
	181	AliError(Form("ITSSA vertex histograms have too few entries for fitting"));
ed15b00c	182	useITSSAvtx=kFALSE;
ed15b00c	183	useSPDvtx=kTRUE;
adca66d2	184	}
adca66d2	185	}
be5ba8df	186
	187	Float_t nEntriesX = histSPDvtxX->GetEffectiveEntries();
	188	Float_t nEntriesY = histSPDvtxY->GetEffectiveEntries();
	189	Float_t nEntriesZ = histSPDvtxZ->GetEffectiveEntries();
e0ea2e34	190
be5ba8df	191	if (nEntriesX < 50. \|\| nEntriesY<50. \|\| nEntriesZ<50.) {
	192	spdAvailable = kFALSE;
	193	if ((useTRKvtx==kFALSE) && (useITSSAvtx==kFALSE)){
	194	AliError(Form("Also SPD vertex histograms have too few entries for fitting, return"));
	195	return;
adca66d2	196	}
be5ba8df	197	}
e0ea2e34	198
a680415a	199	if((nEntriesX == 0.)&&(nEntriesY==0.) && (nEntriesZ>0.)){
	200	vertexerSPD3Doff = kTRUE;
	201	AliWarning("Vertexer SPD 3D off");
	202	}
e0ea2e34	203
	204	Double_t xMeanVtx=0., yMeanVtx=0., zMeanVtx=0.;
	205	Double_t xSigmaVtx=0., ySigmaVtx=0., zSigmaVtx=0.;
	206
be5ba8df	207
e0ea2e34	208	TF1 fitVtxX, fitVtxY, *fitVtxZ;
	209
	210	if (useTRKvtx \|\| useITSSAvtx){
	211	histTRKvtxX ->Fit("gaus", "M");
	212	fitVtxX = histTRKvtxX -> GetFunction("gaus");
	213	xMeanVtx = fitVtxX -> GetParameter(1);
	214	if (TMath::Abs(xMeanVtx) > 2.) {
	215	xMeanVtx = 0.;
be5ba8df	216	writeMeanVertexSPD=kTRUE;
e0ea2e34	217	}
	218
	219	histTRKvtxY ->Fit("gaus", "M");
	220	fitVtxY = histTRKvtxY -> GetFunction("gaus");
	221	yMeanVtx = fitVtxY -> GetParameter(1);
	222	if (TMath::Abs(yMeanVtx) > 2.) {
	223	yMeanVtx = 0.;
be5ba8df	224	writeMeanVertexSPD=kTRUE;
e0ea2e34	225	}
	226
	227	histTRKvtxZ ->Fit("gaus", "M");
	228	fitVtxZ = histTRKvtxZ -> GetFunction("gaus");
	229	zMeanVtx = fitVtxZ -> GetParameter(1);
	230	zSigmaVtx = fitVtxZ -> GetParameter(2);
be5ba8df	231	if ((TMath::Abs(zMeanVtx) > 20.) \|\| (zSigmaVtx>12.)) {
f85d4cb9	232	zMeanVtx = histTRKvtxZ->GetMean();
	233	zSigmaVtx = histTRKvtxZ->GetRMS();
	234	writeMeanVertexSPD=kTRUE;
e0ea2e34	235	}
f85d4cb9	236
e0ea2e34	237	}
	238
	239
	240	//check fits: compare histo mean with fit mean value
	241	Double_t xHistoMean, yHistoMean, zHistoMean;
	242	Double_t xHistoRMS, yHistoRMS, zHistoRMS;
	243
	244	if (useTRKvtx \|\| useITSSAvtx){
f85d4cb9	245	xHistoMean = histTRKvtxX -> GetMean();
	246	xHistoRMS = histTRKvtxX ->GetRMS();
	247
	248	if ((TMath::Abs(xHistoMean-xMeanVtx) > 0.5)){
	249	AliWarning(Form("Possible problems with the fit mean very different from histo mean... using SPD vertex"));
e0ea2e34	250	useTRKvtx = kFALSE;
f85d4cb9	251	useITSSAvtx = kFALSE;
	252	useSPDvtx = kTRUE;
	253	}
	254
	255	yHistoMean = histTRKvtxY ->GetMean();
	256	yHistoRMS = histTRKvtxY ->GetRMS();
	257
	258	if ((TMath::Abs(yHistoMean-yMeanVtx) > 0.5)){
	259	AliWarning(Form("Possible problems with the fit mean very different from histo mean... using SPD vertex"));
e0ea2e34	260	useTRKvtx = kFALSE;
f85d4cb9	261	useITSSAvtx = kFALSE;
	262	useSPDvtx = kTRUE;
	263	}
	264
	265	zHistoMean = histTRKvtxZ -> GetMean();
	266	zHistoRMS = histTRKvtxZ ->GetRMS();
	267
	268	if ((TMath::Abs(zHistoMean-zMeanVtx) > 1.)){
	269	AliWarning(Form("Possible problems with the fit mean very different from histo mean... using SPD vertex"));
e0ea2e34	270	useTRKvtx = kFALSE;
f85d4cb9	271	useITSSAvtx = kFALSE;
	272	useSPDvtx = kTRUE;
	273	}
e0ea2e34	274	}
e0ea2e34	275
f85d4cb9	276
a680415a	277	if ((useSPDvtx) && (spdAvailable) && (!vertexerSPD3Doff)){
f85d4cb9	278
	279	histSPDvtxX ->Fit("gaus", "M");
	280	fitVtxX = histSPDvtxX -> GetFunction("gaus");
	281	xMeanVtx = fitVtxX -> GetParameter(1);
	282	if (TMath::Abs(xMeanVtx) > 2.) {
e0ea2e34	283	xMeanVtx = 0.;
be5ba8df	284	writeMeanVertexSPD=kTRUE;
f85d4cb9	285	}
e0ea2e34	286
f85d4cb9	287	histSPDvtxY ->Fit("gaus", "M");
	288	fitVtxY = histSPDvtxY -> GetFunction("gaus");
	289	yMeanVtx = fitVtxY -> GetParameter(1);
e0ea2e34	290	if (TMath::Abs(yMeanVtx) > 2.) {
e0ea2e34	291	yMeanVtx = 0.;
be5ba8df	292	writeMeanVertexSPD=kTRUE;
e0ea2e34	293	}
	294
	295	histSPDvtxZ ->Fit("gaus", "M");
	296	fitVtxZ = histSPDvtxZ -> GetFunction("gaus");
	297	zMeanVtx = fitVtxZ -> GetParameter(1);
	298	zSigmaVtx = fitVtxZ -> GetParameter(2);
be5ba8df	299	if ((TMath::Abs(zMeanVtx) > 20.) \|\| (zSigmaVtx>12.)) {
f85d4cb9	300	zMeanVtx = histSPDvtxZ ->GetMean();
	301	zSigmaVtx = histSPDvtxZ->GetRMS();
	302	writeMeanVertexSPD = kTRUE;
e0ea2e34	303	}
	304
	305	}
a680415a	306	else if ((useSPDvtx) && (!spdAvailable)){
a680415a	307	AliError(Form("Difference between trkVtx and online one, SPD histos not enough entry or SPD 3D vertex off. Writing Mean Vertex SPD"));
be5ba8df	308	writeMeanVertexSPD = kTRUE;
be5ba8df	309	}
e0ea2e34	310
	311
	312	//check with online position
	313
be5ba8df	314	Double_t posOnline[3], sigmaOnline[3];
be5ba8df	315
f85d4cb9	316	if (useTRKvtx \|\| useITSSAvtx \|\| writeMeanVertexSPD){
e0ea2e34	317	AliCDBManager *manCheck = AliCDBManager::Instance();
e0ea2e34	318	manCheck->SetDefaultStorage("raw://");
adca66d2	319	manCheck->SetRun(runNb);
e0ea2e34	320
e0ea2e34	321	AliCDBEntry *entr = manCheck->Get("GRP/Calib/MeanVertexSPD");
adca66d2	322	if(entr) {
adca66d2	323	AliESDVertex vtxOnline = (AliESDVertex)entr->GetObject();
e0ea2e34	324
adca66d2	325	posOnline[0] = vtxOnline->GetX();
	326	posOnline[1] = vtxOnline->GetY();
	327	posOnline[2] = vtxOnline->GetZ();
e0ea2e34	328
adca66d2	329	sigmaOnline[0] = vtxOnline->GetXRes();
	330	sigmaOnline[1] = vtxOnline->GetYRes();
	331	sigmaOnline[2] = vtxOnline->GetZRes();
e0ea2e34	332
adca66d2	333	//vtxOnline->GetSigmaXYZ(sigmaOnline);
e0ea2e34	334
adca66d2	335	if ((TMath::Abs(posOnline[0]-xMeanVtx) > 0.1) \|\| (TMath::Abs(posOnline[1]-yMeanVtx) > 0.1) \|\| (TMath::Abs(posOnline[2]-zMeanVtx) > 1.)){
adca66d2	336	AliWarning(Form("vertex offline far from the online one"));
adca66d2	337	}
e0ea2e34	338	}
e0ea2e34	339	}
e0ea2e34	340
be5ba8df	341
	342
	343	if (writeMeanVertexSPD){
	344
	345	AliWarning(Form("Writing Mean Vertex SPD, Mean Vertex not available"));
	346
f85d4cb9	347	Double_t sigma[3]={0.0150, 0.0150, zSigmaVtx};
be5ba8df	348
	349	AliESDVertex *vertex = new AliESDVertex(posOnline, sigma, "vertex");
	350
	351	AliCDBManager *cdb = AliCDBManager::Instance();
	352	AliCDBStorage *sto = cdb->GetStorage(dbString);
	353
	354
	355	if (!sto) {
	356	AliError(Form("cannot get storage %s", dbString));
	357	return;
	358	}
	359
	360	AliCDBId id("GRP/Calib/MeanVertex", runNb, runNb);
	361
	362	AliCDBMetaData metaData;
	363	metaData.SetBeamPeriod(0); //check!!!!
	364	metaData.SetResponsible("Davide Caffarri");
	365	metaData.SetComment("Mean Vertex object used in reconstruction");
	366
	367	if (!sto->Put(vertex, id, &metaData)) {
	368	AliError(Form("Error while putting object in storage %s", dbString));
	369	}
	370
	371	delete vertex;
	372	return;
	373	}
	374
	375
e0ea2e34	376	Float_t meanMult = 40.;
	377	Float_t p2 = 1.4;
	378	Float_t resolVtx = 0.05;
	379
	380	Double_t xSigmaMult, ySigmaMult, corrXZ, corrYZ, lumiRegSquaredX, lumiRegSquaredY;
	381	Double_t covarXZ=0., covarYZ=0.;
e0ea2e34	382
a680415a	383	Bool_t highMultEnvironment = kTRUE;
a680415a	384
e0ea2e34	385	TF1 sigmaFitX, sigmaFitY, *corrFit;
	386
	387	TH1F *histTRKdefMultX=0;
	388	TH1F *histTRKdefMultY=0;
a680415a	389	TH1F *histTRKHighMultX=0;
a680415a	390	TH1F *histTRKHighMultY=0;
e0ea2e34	391	TH2F *histTRKVertexXZ=0;
	392	TH2F *histTRKVertexYZ=0;
	393
e0ea2e34	394
adca66d2	395	if (useTRKvtx){
792b5897	396	if (list){
	397	histTRKdefMultX = (TH1F*)list->FindObject("hTRKVertexXdefMult");
	398	histTRKdefMultY = (TH1F*)list->FindObject("hTRKVertexYdefMult");
a680415a	399	histTRKHighMultX = (TH1F*)list->FindObject("hTRKVertexXHighMult");
a680415a	400	histTRKHighMultY = (TH1F*)list->FindObject("hTRKVertexYHighMult");
792b5897	401	histTRKVertexXZ = (TH2F*)list->FindObject("hTRKVertexXZ");
	402	histTRKVertexYZ = (TH2F*)list->FindObject("hTRKVertexYZ");
	403	}
	404
	405	else {
	406	histTRKdefMultX = (TH1F*)file->Get("hTRKVertexXdefMult");
	407	histTRKdefMultY = (TH1F*)file->Get("hTRKVertexYdefMult");
a680415a	408	histTRKHighMultX = (TH1F*)file->Get("hTRKVertexXHighMult");
a680415a	409	histTRKHighMultY = (TH1F*)file->Get("hTRKVertexYHighMult");
792b5897	410	histTRKVertexXZ = (TH2F*)file->Get("hTRKVertexXZ");
	411	histTRKVertexYZ = (TH2F*)file->Get("hTRKVertexYZ");
	412	}
	413
adca66d2	414	}
e0ea2e34	415
adca66d2	416	if (useITSSAvtx){
792b5897	417	if (list){
	418	histTRKdefMultX = (TH1F*)list->FindObject("hITSSAVertexXdefMult");
	419	histTRKdefMultY = (TH1F*)list->FindObject("hITSSAVertexYdefMult");
a680415a	420	histTRKHighMultX = (TH1F*)list->FindObject("hITSSAVertexXHighMult");
a680415a	421	histTRKHighMultY = (TH1F*)list->FindObject("hITSSAVertexYHighMult");
792b5897	422	histTRKVertexXZ = (TH2F*)list->FindObject("hITSSAVertexXZ");
	423	histTRKVertexYZ = (TH2F*)list->FindObject("hITSSAVertexYZ");
	424	}
adca66d2	425
792b5897	426	else {
	427	histTRKdefMultX = (TH1F*)file->Get("hITSSAVertexXdefMult");
	428	histTRKdefMultY = (TH1F*)file->Get("hITSSAVertexYdefMult");
a680415a	429	histTRKHighMultX = (TH1F*)file->Get("hITSSAVertexXHighMult");
a680415a	430	histTRKHighMultY = (TH1F*)file->Get("hITSSAVertexYHighMult");
792b5897	431	histTRKVertexXZ = (TH2F*)file->Get("hITSSAVertexXZ");
	432	histTRKVertexYZ = (TH2F*)file->Get("hITSSAVertexYZ");
	433	}
a680415a	434	}
792b5897	435
a680415a	436
a680415a	437	if ((histTRKHighMultX) && (histTRKHighMultY)){
792b5897	438
a680415a	439	Float_t nEntriesHighMultX = histTRKHighMultX->GetEffectiveEntries();
a680415a	440	Float_t nEntriesHighMultY = histTRKHighMultY->GetEffectiveEntries();
e0ea2e34	441
a680415a	442	if ((nEntriesHighMultX >100) && (nEntriesHighMultY>100)) {
	443	AliWarning(Form("Setting High Mulitplicity environment"));
	444	highMultEnvironment = kTRUE;
	445	}
	446	else {
	447	AliWarning(Form("Setting Low Mulitplicity environment"));
	448	highMultEnvironment=kFALSE;
	449	}
	450	}
e0ea2e34	451
a680415a	452	else{
	453	AliWarning(Form("No histos found, setting Low Mulitplicity environment"));
	454	highMultEnvironment=kFALSE;
	455
	456	}
e0ea2e34	457
a680415a	458	if (highMultEnvironment==kFALSE){
adca66d2	459
a680415a	460	if ((!histTRKdefMultX) \|\| (!histTRKdefMultY) \|\| (histTRKdefMultX->GetEntries() < 40.) \|\| (histTRKdefMultY->GetEntries() < 40.)){
	461	AliWarning(Form("histos for lumi reg calculation not found, default value setted"));
	462	xSigmaVtx=0.0120;
	463	ySigmaVtx=0.0120;
	464	} else {
	465
	466	histTRKdefMultX -> Fit("gaus", "M", "", -0.4, 0.4);
	467	sigmaFitX = histTRKdefMultX -> GetFunction("gaus");
	468	xSigmaMult = sigmaFitX->GetParameter(2);
adca66d2	469
a680415a	470	lumiRegSquaredX = (xSigmaMultxSigmaMult - ((resolVtxresolVtx)/TMath::Power(meanMult, p2)));
	471	if (lumiRegSquaredX < 0) {
	472	AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
	473	xSigmaMult = 0.;
	474	xSigmaVtx = 0.0120;
	475	}
adca66d2	476
a680415a	477	if (lumiRegSquaredX > 0 && lumiRegSquaredX < 0.0005){
	478	xSigmaVtx = TMath::Sqrt(lumiRegSquaredX);
	479	xSigmaVtx = xSigmaVtx*1.2;
	480	}
adca66d2	481
a680415a	482	histTRKdefMultY -> Fit("gaus", "M", "", -0.2, 0.6);
	483	sigmaFitY = histTRKdefMultY -> GetFunction("gaus");
	484	ySigmaMult = sigmaFitY->GetParameter(2);
adca66d2	485
a680415a	486	lumiRegSquaredY= (ySigmaMultySigmaMult - ((resolVtxresolVtx)/TMath::Power(meanMult, p2)));
	487	if (lumiRegSquaredY < 0) {
	488	AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
	489	ySigmaMult = 0.;
	490	ySigmaVtx = 0.0120;
	491	}
adca66d2	492
a680415a	493	if (lumiRegSquaredY > 0 && lumiRegSquaredY < 0.0005){
	494	ySigmaVtx = TMath::Sqrt(lumiRegSquaredY);
	495	ySigmaVtx = ySigmaVtx*1.2;
	496	}
	497
	498	TProfile *htrkXZ = histTRKVertexXZ ->ProfileY();
	499	htrkXZ -> Fit("pol1", "M", "", -10., 10.);
	500	corrFit = htrkXZ->GetFunction("pol1");
	501	corrXZ = corrFit->GetParameter(1);
	502
	503	if (TMath::Abs(corrXZ) > 0.01) {
	504	AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
	505	corrXZ =0.;
	506	}
	507	else{
	508	covarXZ = corrXZ * zSigmaVtx*zSigmaVtx;
	509
	510	}
	511
	512	TProfile *htrkYZ = histTRKVertexYZ ->ProfileY();
	513	htrkYZ -> Fit("pol1", "M", "", -10., 10.);
	514	corrFit = htrkYZ->GetFunction("pol1");
	515	corrYZ = corrFit->GetParameter(1);
	516
	517	if (TMath::Abs(corrYZ) > 0.01) {
	518	AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
	519	corrYZ =0.;
	520	}
	521	else{
	522	covarYZ = corrYZzSigmaVtxzSigmaVtx;
	523	}
	524
	525	}
	526	}
adca66d2	527
a680415a	528	if (highMultEnvironment==kTRUE){
	529
	530	histTRKHighMultX -> Fit("gaus", "M", "", -0.4, 0.4);
	531	sigmaFitX = histTRKHighMultX -> GetFunction("gaus");
	532	xSigmaMult = sigmaFitX->GetParameter(2);
	533
	534	if ((xSigmaMult <0) \|\| (xSigmaMult>0.03)){
	535	AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
	536	xSigmaMult = 0.;
	537	xSigmaVtx = 0.0120;
	538	}
	539	else{
	540	xSigmaVtx = xSigmaMult;
	541	xSigmaVtx = xSigmaVtx*1.2;
	542	}
	543
	544	histTRKHighMultY -> Fit("gaus", "M", "", -0.2, 0.5);
	545	sigmaFitY = histTRKHighMultY -> GetFunction("gaus");
	546	ySigmaMult = sigmaFitY->GetParameter(2);
	547
	548	if ((ySigmaMult <0) \|\| (ySigmaMult>0.03)){
	549	AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
	550	ySigmaMult = 0.;
	551	ySigmaVtx = 0.0120;
	552	}
	553	else{
	554	ySigmaVtx = ySigmaMult;
	555	ySigmaVtx = ySigmaVtx*1.2;
	556	}
	557
	558	TProfile *htrkXZ = histTRKVertexXZ ->ProfileY();
	559	htrkXZ -> Fit("pol1", "M", "", -10., 10.);
	560	corrFit = htrkXZ->GetFunction("pol1");
	561	corrXZ = corrFit->GetParameter(1);
adca66d2	562
	563	if (TMath::Abs(corrXZ) > 0.01) {
	564	AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
	565	corrXZ =0.;
	566	}
	567	else{
	568	covarXZ = corrXZ * zSigmaVtx*zSigmaVtx;
	569
	570	}
a680415a	571
adca66d2	572	TProfile *htrkYZ = histTRKVertexYZ ->ProfileY();
	573	htrkYZ -> Fit("pol1", "M", "", -10., 10.);
	574	corrFit = htrkYZ->GetFunction("pol1");
	575	corrYZ = corrFit->GetParameter(1);
	576
	577	if (TMath::Abs(corrYZ) > 0.01) {
	578	AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
	579	corrYZ =0.;
	580	}
	581	else{
	582	covarYZ = corrYZzSigmaVtxzSigmaVtx;
	583	}
	584
e0ea2e34	585	}
a680415a	586
e0ea2e34	587	Double_t position[3], covMatrix[6];
	588	Double_t chi2=1.;
	589	Int_t nContr=1;
	590
	591	position[0] = xMeanVtx;
	592	position[1] = yMeanVtx;
	593	position[2] = zMeanVtx;
	594
792b5897	595	covMatrix[0] = xSigmaVtx*xSigmaVtx;
e0ea2e34	596	covMatrix[1] = 0.; //xy
792b5897	597	covMatrix[2] = ySigmaVtx*ySigmaVtx;
e0ea2e34	598	covMatrix[3] = covarXZ;
e0ea2e34	599	covMatrix[4] = covarYZ;
792b5897	600	covMatrix[5] = zSigmaVtx*zSigmaVtx;
e0ea2e34	601
e0ea2e34	602
792b5897	603	//Printf ("%f, %f, %f, %f", xSigmaVtx, ySigmaVtx, covarXZ, covarYZ);
e0ea2e34	604
	605	AliESDVertex *vertex = new AliESDVertex(position, covMatrix, chi2, nContr, "vertex");
	606
792b5897	607	AliCDBManager *cdb = AliCDBManager::Instance();
e0ea2e34	608	AliCDBStorage *sto = cdb->GetStorage(dbString);
792b5897	609
	610
	611	if (!sto) {
e0ea2e34	612	AliError(Form("cannot get storage %s", dbString));
	613	return;
	614	}
	615
	616	AliCDBId id("GRP/Calib/MeanVertex", runNb, runNb);
	617
	618	AliCDBMetaData metaData;
	619	metaData.SetBeamPeriod(0); //check!!!!
	620	metaData.SetResponsible("Davide Caffarri");
	621	metaData.SetComment("Mean Vertex object used in reconstruction");
	622
792b5897	623	if (!sto->Put(vertex, id, &metaData)) {
792b5897	624	AliError(Form("Error while putting object in storage %s", dbString));
e0ea2e34	625	}
	626
	627	delete vertex;
	628
	629
	630
	631	}
	632
	633