PWGPP/ITS/AliMeanVertexPreprocessorOffline.cxx

   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
  71     if (!dbString){
  72                 AliError(Form("no OCDB path found, return"));
  73             return;
  74         }
  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
  87
  88         Bool_t useTRKvtx = kTRUE;
  89         Bool_t useITSSAvtx = kFALSE;
  90         Bool_t useSPDvtx = kFALSE;
  91         Bool_t spdAvailable = kTRUE;
  92         Bool_t writeMeanVertexSPD = kFALSE;
  93         Bool_t vertexerSPD3Doff=kFALSE;
  94
  95
  96     if (!list) {
  97
  98                 histTRKvtxX = (TH1F*)file->Get("hTRKVertexX");
  99                 histTRKvtxY = (TH1F*)file->Get("hTRKVertexY");
 100                 histTRKvtxZ = (TH1F*)file->Get("hTRKVertexZ");
 101
 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");
 114
 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{
 129
 130                 histTRKvtxX = (TH1F*)list->FindObject("hTRKVertexX");
 131                 histTRKvtxY = (TH1F*)list->FindObject("hTRKVertexY");
 132                 histTRKvtxZ = (TH1F*)list->FindObject("hTRKVertexZ");
 133
 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;
 142
 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         }
 160
 161
 162         if (useTRKvtx){
 163
 164                 Float_t nEntriesX = histTRKvtxX->GetEffectiveEntries();
 165                 Float_t nEntriesY = histTRKvtxY->GetEffectiveEntries();
 166                 Float_t nEntriesZ = histTRKvtxZ->GetEffectiveEntries();
 167
 168                 if (nEntriesX < 50. || nEntriesY<50. || nEntriesZ<50.) {
 169                         AliError(Form("TRK vertex histograms have too few entries for fitting"));
 170                         useTRKvtx=kFALSE;
 171                         useSPDvtx = kTRUE;
 172                 }
 173         }
 174         if (useITSSAvtx){
 175
 176                 Float_t nEntriesX = histTRKvtxX->GetEffectiveEntries();
 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"));
 182                         useITSSAvtx=kFALSE;
 183                         useSPDvtx=kTRUE;
 184                 }
 185         }
 186
 187         Float_t nEntriesX = histSPDvtxX->GetEffectiveEntries();
 188         Float_t nEntriesY = histSPDvtxY->GetEffectiveEntries();
 189         Float_t nEntriesZ = histSPDvtxZ->GetEffectiveEntries();
 190
 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;
 196                 }
 197         }
 198
 199         if((nEntriesX == 0.)&&(nEntriesY==0.) && (nEntriesZ>0.)){
 200                 vertexerSPD3Doff = kTRUE;
 201                 AliWarning("Vertexer SPD 3D off");
 202         }
 203
 204         Double_t xMeanVtx=0., yMeanVtx=0., zMeanVtx=0.;
 205         Double_t xSigmaVtx=0., ySigmaVtx=0., zSigmaVtx=0.;
 206
 207
 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.;
 216                         writeMeanVertexSPD=kTRUE;
 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.;
 224                         writeMeanVertexSPD=kTRUE;
 225                 }
 226
 227                 histTRKvtxZ ->Fit("gaus", "M");
 228                 fitVtxZ = histTRKvtxZ -> GetFunction("gaus");
 229                 zMeanVtx = fitVtxZ -> GetParameter(1);
 230                 zSigmaVtx = fitVtxZ -> GetParameter(2);
 231                 if ((TMath::Abs(zMeanVtx) > 20.) || (zSigmaVtx>12.)) {
 232                         zMeanVtx = 0.;
 233                         zMeanVtx = 5.;
 234                         writeMeanVertexSPD=kTRUE;
 235                 }
 236
 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){
 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"));
 250             useTRKvtx = kFALSE;
 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"));
 260             useTRKvtx = kFALSE;
 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"));
 270             useTRKvtx = kFALSE;
 271                         useITSSAvtx = kFALSE;
 272                         useSPDvtx = kTRUE;
 273                 }
 274         }
 275
 276
 277         if ((useSPDvtx) && (spdAvailable) && (!vertexerSPD3Doff)){
 278
 279                 histSPDvtxX ->Fit("gaus", "M");
 280                 fitVtxX = histSPDvtxX -> GetFunction("gaus");
 281                 xMeanVtx = fitVtxX -> GetParameter(1);
 282                 if (TMath::Abs(xMeanVtx) > 2.) {
 283                         xMeanVtx = 0.;
 284                         writeMeanVertexSPD=kTRUE;
 285                 }
 286
 287                 histSPDvtxY ->Fit("gaus", "M");
 288                 fitVtxY = histSPDvtxY -> GetFunction("gaus");
 289                 yMeanVtx = fitVtxY -> GetParameter(1);
 290                 if (TMath::Abs(yMeanVtx) > 2.) {
 291                         yMeanVtx = 0.;
 292                         writeMeanVertexSPD=kTRUE;
 293                 }
 294
 295                 histSPDvtxZ ->Fit("gaus", "M");
 296                 fitVtxZ = histSPDvtxZ -> GetFunction("gaus");
 297                 zMeanVtx = fitVtxZ -> GetParameter(1);
 298                 zSigmaVtx = fitVtxZ -> GetParameter(2);
 299                 if ((TMath::Abs(zMeanVtx) > 20.) || (zSigmaVtx>12.)) {
 300                         zMeanVtx = 0.;
 301                         zSigmaVtx = 5.;
 302                         writeMeanVertexSPD = kTRUE;
 303                 }
 304
 305         }
 306         else if ((useSPDvtx) && (!spdAvailable)){
 307                 AliError(Form("Difference between trkVtx and online one, SPD histos not enough entry or SPD 3D vertex off. Writing Mean Vertex SPD"));
 308                 writeMeanVertexSPD = kTRUE;
 309         }
 310
 311
 312         //check with online position
 313
 314          Double_t posOnline[3], sigmaOnline[3];
 315
 316         if (useTRKvtx || useITSSAvtx){
 317                 AliCDBManager *manCheck = AliCDBManager::Instance();
 318                 manCheck->SetDefaultStorage("raw://");
 319                 manCheck->SetRun(runNb);
 320
 321                 AliCDBEntry *entr = manCheck->Get("GRP/Calib/MeanVertexSPD");
 322                 if(entr) {
 323                   AliESDVertex *vtxOnline = (AliESDVertex*)entr->GetObject();
 324
 325                   posOnline[0] = vtxOnline->GetX();
 326                   posOnline[1] = vtxOnline->GetY();
 327                   posOnline[2] = vtxOnline->GetZ();
 328
 329                   sigmaOnline[0] = vtxOnline->GetXRes();
 330                   sigmaOnline[1] = vtxOnline->GetYRes();
 331                   sigmaOnline[2] = vtxOnline->GetZRes();
 332
 333                   //vtxOnline->GetSigmaXYZ(sigmaOnline);
 334
 335                   if ((TMath::Abs(posOnline[0]-xMeanVtx) > 0.1) || (TMath::Abs(posOnline[1]-yMeanVtx) > 0.1) || (TMath::Abs(posOnline[2]-zMeanVtx) > 1.)){
 336                     AliWarning(Form("vertex offline far from the online one"));
 337                   }
 338                 }
 339         }
 340
 341
 342
 343         if (writeMeanVertexSPD){
 344
 345                 AliWarning(Form("Writing Mean Vertex SPD, Mean Vertex not available"));
 346
 347                 Double_t sigma[3]={0.0150, 0.0150, 6.};
 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
 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.;
 382
 383         Bool_t highMultEnvironment = kTRUE;
 384
 385         TF1 *sigmaFitX, *sigmaFitY, *corrFit;
 386
 387         TH1F *histTRKdefMultX=0;
 388         TH1F *histTRKdefMultY=0;
 389         TH1F *histTRKHighMultX=0;
 390         TH1F *histTRKHighMultY=0;
 391         TH2F *histTRKVertexXZ=0;
 392         TH2F *histTRKVertexYZ=0;
 393
 394
 395         if (useTRKvtx){
 396           if (list){
 397                   histTRKdefMultX = (TH1F*)list->FindObject("hTRKVertexXdefMult");
 398                   histTRKdefMultY = (TH1F*)list->FindObject("hTRKVertexYdefMult");
 399                   histTRKHighMultX = (TH1F*)list->FindObject("hTRKVertexXHighMult");
 400                   histTRKHighMultY = (TH1F*)list->FindObject("hTRKVertexYHighMult");
 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");
 408                   histTRKHighMultX = (TH1F*)file->Get("hTRKVertexXHighMult");
 409                   histTRKHighMultY = (TH1F*)file->Get("hTRKVertexYHighMult");
 410                   histTRKVertexXZ = (TH2F*)file->Get("hTRKVertexXZ");
 411                   histTRKVertexYZ = (TH2F*)file->Get("hTRKVertexYZ");
 412           }
 413
 414         }
 415
 416         if (useITSSAvtx){
 417           if (list){
 418                   histTRKdefMultX = (TH1F*)list->FindObject("hITSSAVertexXdefMult");
 419                   histTRKdefMultY = (TH1F*)list->FindObject("hITSSAVertexYdefMult");
 420                   histTRKHighMultX = (TH1F*)list->FindObject("hITSSAVertexXHighMult");
 421                   histTRKHighMultY = (TH1F*)list->FindObject("hITSSAVertexYHighMult");
 422                   histTRKVertexXZ = (TH2F*)list->FindObject("hITSSAVertexXZ");
 423                   histTRKVertexYZ = (TH2F*)list->FindObject("hITSSAVertexYZ");
 424           }
 425
 426           else {
 427                   histTRKdefMultX = (TH1F*)file->Get("hITSSAVertexXdefMult");
 428                   histTRKdefMultY = (TH1F*)file->Get("hITSSAVertexYdefMult");
 429                   histTRKHighMultX = (TH1F*)file->Get("hITSSAVertexXHighMult");
 430                   histTRKHighMultY = (TH1F*)file->Get("hITSSAVertexYHighMult");
 431                   histTRKVertexXZ = (TH2F*)file->Get("hITSSAVertexXZ");
 432                   histTRKVertexYZ = (TH2F*)file->Get("hITSSAVertexYZ");
 433           }
 434         }
 435
 436
 437         if ((histTRKHighMultX) && (histTRKHighMultY)){
 438
 439                 Float_t nEntriesHighMultX = histTRKHighMultX->GetEffectiveEntries();
 440                 Float_t nEntriesHighMultY = histTRKHighMultY->GetEffectiveEntries();
 441
 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         }
 451
 452         else{
 453                 AliWarning(Form("No histos found, setting Low Mulitplicity environment"));
 454                 highMultEnvironment=kFALSE;
 455
 456         }
 457
 458         if (highMultEnvironment==kFALSE){
 459
 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);
 469
 470                         lumiRegSquaredX = (xSigmaMult*xSigmaMult - ((resolVtx*resolVtx)/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                         }
 476
 477                         if (lumiRegSquaredX > 0 && lumiRegSquaredX < 0.0005){
 478                                 xSigmaVtx = TMath::Sqrt(lumiRegSquaredX);
 479                                 xSigmaVtx = xSigmaVtx*1.2;
 480                         }
 481
 482                         histTRKdefMultY -> Fit("gaus", "M", "", -0.2, 0.6);
 483                         sigmaFitY = histTRKdefMultY -> GetFunction("gaus");
 484                         ySigmaMult = sigmaFitY->GetParameter(2);
 485
 486                         lumiRegSquaredY= (ySigmaMult*ySigmaMult - ((resolVtx*resolVtx)/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                         }
 492
 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 = corrYZ*zSigmaVtx*zSigmaVtx;
 523                         }
 524
 525                 }
 526         }
 527
 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);
 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           }
 571
 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 = corrYZ*zSigmaVtx*zSigmaVtx;
 583           }
 584
 585         }
 586
 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
 595   covMatrix[0] = xSigmaVtx*xSigmaVtx;
 596   covMatrix[1] = 0.; //xy
 597   covMatrix[2] = ySigmaVtx*ySigmaVtx;
 598   covMatrix[3] = covarXZ;
 599   covMatrix[4] = covarYZ;
 600   covMatrix[5] = zSigmaVtx*zSigmaVtx;
 601
 602
 603         //Printf ("%f, %f, %f, %f", xSigmaVtx, ySigmaVtx, covarXZ, covarYZ);
 604
 605   AliESDVertex  *vertex =  new AliESDVertex(position, covMatrix, chi2, nContr, "vertex");
 606
 607   AliCDBManager *cdb = AliCDBManager::Instance();
 608   AliCDBStorage *sto = cdb->GetStorage(dbString);
 609
 610
 611         if (!sto) {
 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
 623   if (!sto->Put(vertex, id, &metaData)) {
 624         AliError(Form("Error while putting object in storage %s", dbString));
 625    }
 626
 627   delete vertex;
 628
 629
 630
 631 }
 632
 633