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, AliCDBStorage *db, 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 (!db){
  72           AliError(Form("no OCDB storage 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 = histTRKvtxZ->GetMean();
 233                   zSigmaVtx = histTRKvtxZ->GetRMS();
 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 = histSPDvtxZ ->GetMean();
 301                   zSigmaVtx = histSPDvtxZ->GetRMS();
 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 || writeMeanVertexSPD){
 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, zSigmaVtx};
 348
 349                 AliESDVertex  *vertex =  new AliESDVertex(posOnline, sigma, "vertex");
 350
 351                 AliCDBId id("GRP/Calib/MeanVertex", runNb, runNb);
 352
 353                 AliCDBMetaData metaData;
 354                 metaData.SetBeamPeriod(0); //check!!!!
 355                 metaData.SetResponsible("Davide Caffarri");
 356                 metaData.SetComment("Mean Vertex object used in reconstruction");
 357
 358                 if (!db->Put(vertex, id, &metaData)) {
 359                   AliError(Form("Error while putting object in storage %s", db->GetURI().Data()));
 360                 }
 361
 362                 delete vertex;
 363                 return;
 364         }
 365
 366
 367         Float_t meanMult = 40.;
 368         Float_t p2 = 1.4;
 369         Float_t resolVtx = 0.05;
 370
 371         Double_t xSigmaMult, ySigmaMult, corrXZ, corrYZ, lumiRegSquaredX, lumiRegSquaredY;
 372         Double_t covarXZ=0., covarYZ=0.;
 373
 374         Bool_t highMultEnvironment = kTRUE;
 375
 376         TF1 *sigmaFitX, *sigmaFitY, *corrFit;
 377
 378         TH1F *histTRKdefMultX=0;
 379         TH1F *histTRKdefMultY=0;
 380         TH1F *histTRKHighMultX=0;
 381         TH1F *histTRKHighMultY=0;
 382         TH2F *histTRKVertexXZ=0;
 383         TH2F *histTRKVertexYZ=0;
 384
 385
 386         if (useTRKvtx){
 387           if (list){
 388                   histTRKdefMultX = (TH1F*)list->FindObject("hTRKVertexXdefMult");
 389                   histTRKdefMultY = (TH1F*)list->FindObject("hTRKVertexYdefMult");
 390                   histTRKHighMultX = (TH1F*)list->FindObject("hTRKVertexXHighMult");
 391                   histTRKHighMultY = (TH1F*)list->FindObject("hTRKVertexYHighMult");
 392                   histTRKVertexXZ = (TH2F*)list->FindObject("hTRKVertexXZ");
 393                   histTRKVertexYZ = (TH2F*)list->FindObject("hTRKVertexYZ");
 394           }
 395
 396           else {
 397                   histTRKdefMultX = (TH1F*)file->Get("hTRKVertexXdefMult");
 398                   histTRKdefMultY = (TH1F*)file->Get("hTRKVertexYdefMult");
 399                   histTRKHighMultX = (TH1F*)file->Get("hTRKVertexXHighMult");
 400                   histTRKHighMultY = (TH1F*)file->Get("hTRKVertexYHighMult");
 401                   histTRKVertexXZ = (TH2F*)file->Get("hTRKVertexXZ");
 402                   histTRKVertexYZ = (TH2F*)file->Get("hTRKVertexYZ");
 403           }
 404
 405         }
 406
 407         if (useITSSAvtx){
 408           if (list){
 409                   histTRKdefMultX = (TH1F*)list->FindObject("hITSSAVertexXdefMult");
 410                   histTRKdefMultY = (TH1F*)list->FindObject("hITSSAVertexYdefMult");
 411                   histTRKHighMultX = (TH1F*)list->FindObject("hITSSAVertexXHighMult");
 412                   histTRKHighMultY = (TH1F*)list->FindObject("hITSSAVertexYHighMult");
 413                   histTRKVertexXZ = (TH2F*)list->FindObject("hITSSAVertexXZ");
 414                   histTRKVertexYZ = (TH2F*)list->FindObject("hITSSAVertexYZ");
 415           }
 416
 417           else {
 418                   histTRKdefMultX = (TH1F*)file->Get("hITSSAVertexXdefMult");
 419                   histTRKdefMultY = (TH1F*)file->Get("hITSSAVertexYdefMult");
 420                   histTRKHighMultX = (TH1F*)file->Get("hITSSAVertexXHighMult");
 421                   histTRKHighMultY = (TH1F*)file->Get("hITSSAVertexYHighMult");
 422                   histTRKVertexXZ = (TH2F*)file->Get("hITSSAVertexXZ");
 423                   histTRKVertexYZ = (TH2F*)file->Get("hITSSAVertexYZ");
 424           }
 425         }
 426
 427
 428         if ((histTRKHighMultX) && (histTRKHighMultY)){
 429
 430                 Float_t nEntriesHighMultX = histTRKHighMultX->GetEffectiveEntries();
 431                 Float_t nEntriesHighMultY = histTRKHighMultY->GetEffectiveEntries();
 432
 433                 if ((nEntriesHighMultX >100) && (nEntriesHighMultY>100)) {
 434                         AliWarning(Form("Setting High Mulitplicity environment"));
 435                 highMultEnvironment = kTRUE;
 436                 }
 437                 else {
 438                         AliWarning(Form("Setting Low Mulitplicity environment"));
 439                         highMultEnvironment=kFALSE;
 440                 }
 441         }
 442
 443         else{
 444                 AliWarning(Form("No histos found, setting Low Mulitplicity environment"));
 445                 highMultEnvironment=kFALSE;
 446
 447         }
 448
 449         if (highMultEnvironment==kFALSE){
 450
 451                 if ((!histTRKdefMultX) || (!histTRKdefMultY) || (histTRKdefMultX->GetEntries() < 40.) || (histTRKdefMultY->GetEntries() < 40.)){
 452                         AliWarning(Form("histos for lumi reg calculation not found, default value setted"));
 453                         xSigmaVtx=0.0120;
 454                         ySigmaVtx=0.0120;
 455                 } else {
 456
 457                         histTRKdefMultX -> Fit("gaus", "M", "", -0.4, 0.4);
 458                         sigmaFitX = histTRKdefMultX -> GetFunction("gaus");
 459                         xSigmaMult = sigmaFitX->GetParameter(2);
 460
 461                         lumiRegSquaredX = (xSigmaMult*xSigmaMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
 462                         if (lumiRegSquaredX < 0) {
 463                                 AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
 464                                 xSigmaMult = 0.;
 465                                 xSigmaVtx = 0.0120;
 466                         }
 467
 468                         if (lumiRegSquaredX > 0 && lumiRegSquaredX < 0.0005){
 469                                 xSigmaVtx = TMath::Sqrt(lumiRegSquaredX);
 470                                 xSigmaVtx = xSigmaVtx*1.2;
 471                         }
 472
 473                         histTRKdefMultY -> Fit("gaus", "M", "", -0.2, 0.6);
 474                         sigmaFitY = histTRKdefMultY -> GetFunction("gaus");
 475                         ySigmaMult = sigmaFitY->GetParameter(2);
 476
 477                         lumiRegSquaredY= (ySigmaMult*ySigmaMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
 478                         if (lumiRegSquaredY < 0) {
 479                                 AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
 480                                 ySigmaMult = 0.;
 481                                 ySigmaVtx = 0.0120;
 482                         }
 483
 484                         if (lumiRegSquaredY > 0 && lumiRegSquaredY < 0.0005){
 485                                 ySigmaVtx = TMath::Sqrt(lumiRegSquaredY);
 486                                 ySigmaVtx = ySigmaVtx*1.2;
 487                         }
 488
 489                         TProfile *htrkXZ = histTRKVertexXZ ->ProfileY();
 490                         htrkXZ -> Fit("pol1", "M", "", -10., 10.);
 491                         corrFit = htrkXZ->GetFunction("pol1");
 492                         corrXZ = corrFit->GetParameter(1);
 493
 494                         if (TMath::Abs(corrXZ) > 0.01) {
 495                                 AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 496                                 corrXZ =0.;
 497                         }
 498                         else{
 499                                 covarXZ = corrXZ * zSigmaVtx*zSigmaVtx;
 500
 501                         }
 502
 503                         TProfile *htrkYZ = histTRKVertexYZ ->ProfileY();
 504                         htrkYZ -> Fit("pol1", "M", "", -10., 10.);
 505                         corrFit = htrkYZ->GetFunction("pol1");
 506                         corrYZ = corrFit->GetParameter(1);
 507
 508                         if (TMath::Abs(corrYZ) > 0.01) {
 509                                 AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 510                                 corrYZ =0.;
 511                         }
 512                         else{
 513                                 covarYZ = corrYZ*zSigmaVtx*zSigmaVtx;
 514                         }
 515
 516                 }
 517         }
 518
 519         if (highMultEnvironment==kTRUE){
 520
 521                 histTRKHighMultX -> Fit("gaus", "M", "", -0.4, 0.4);
 522                 sigmaFitX = histTRKHighMultX -> GetFunction("gaus");
 523                 xSigmaMult = sigmaFitX->GetParameter(2);
 524
 525                 if ((xSigmaMult <0) || (xSigmaMult>0.03)){
 526                         AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
 527                         xSigmaMult = 0.;
 528                         xSigmaVtx = 0.0120;
 529                 }
 530         else{
 531                         xSigmaVtx = xSigmaMult;
 532                     xSigmaVtx = xSigmaVtx*1.2;
 533                 }
 534
 535                 histTRKHighMultY -> Fit("gaus", "M", "", -0.2, 0.5);
 536                 sigmaFitY = histTRKHighMultY -> GetFunction("gaus");
 537                 ySigmaMult = sigmaFitY->GetParameter(2);
 538
 539                 if ((ySigmaMult <0) || (ySigmaMult>0.03)){
 540                         AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
 541                         ySigmaMult = 0.;
 542                         ySigmaVtx = 0.0120;
 543                 }
 544         else{
 545                         ySigmaVtx = ySigmaMult;
 546                     ySigmaVtx = ySigmaVtx*1.2;
 547                 }
 548
 549            TProfile *htrkXZ = histTRKVertexXZ ->ProfileY();
 550            htrkXZ -> Fit("pol1", "M", "", -10., 10.);
 551            corrFit = htrkXZ->GetFunction("pol1");
 552            corrXZ = corrFit->GetParameter(1);
 553
 554           if (TMath::Abs(corrXZ) > 0.01) {
 555             AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 556             corrXZ =0.;
 557           }
 558           else{
 559             covarXZ = corrXZ * zSigmaVtx*zSigmaVtx;
 560
 561           }
 562
 563           TProfile *htrkYZ = histTRKVertexYZ ->ProfileY();
 564           htrkYZ -> Fit("pol1", "M", "", -10., 10.);
 565           corrFit = htrkYZ->GetFunction("pol1");
 566           corrYZ = corrFit->GetParameter(1);
 567
 568           if (TMath::Abs(corrYZ) > 0.01) {
 569             AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 570             corrYZ =0.;
 571           }
 572           else{
 573             covarYZ = corrYZ*zSigmaVtx*zSigmaVtx;
 574           }
 575
 576         }
 577
 578   Double_t position[3], covMatrix[6];
 579   Double_t chi2=1.;
 580   Int_t nContr=1;
 581
 582   position[0] = xMeanVtx;
 583   position[1] = yMeanVtx;
 584   position[2] = zMeanVtx;
 585
 586   covMatrix[0] = xSigmaVtx*xSigmaVtx;
 587   covMatrix[1] = 0.; //xy
 588   covMatrix[2] = ySigmaVtx*ySigmaVtx;
 589   covMatrix[3] = covarXZ;
 590   covMatrix[4] = covarYZ;
 591   covMatrix[5] = zSigmaVtx*zSigmaVtx;
 592
 593
 594         //Printf ("%f, %f, %f, %f", xSigmaVtx, ySigmaVtx, covarXZ, covarYZ);
 595
 596   AliESDVertex  *vertex =  new AliESDVertex(position, covMatrix, chi2, nContr, "vertex");
 597
 598   AliCDBId id("GRP/Calib/MeanVertex", runNb, runNb);
 599
 600   AliCDBMetaData metaData;
 601   metaData.SetBeamPeriod(0); //check!!!!
 602   metaData.SetResponsible("Davide Caffarri");
 603   metaData.SetComment("Mean Vertex object used in reconstruction");
 604
 605   if (!db->Put(vertex, id, &metaData)) {
 606     AliError(Form("Error while putting object in storage %s", db->GetURI().Data()));
 607    }
 608
 609   delete vertex;
 610
 611
 612
 613 }
 614
 615