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 #include <TCanvas.h>
  39
  40 #include "AliESDVertex.h"
  41 #include "TH1F.h"
  42 #include "TH2F.h"
  43 #include "TF1.h"
  44 #include "TProfile.h"
  45
  46
  47 ClassImp(AliMeanVertexPreprocessorOffline)
  48
  49 //_______________________________________________________
  50
  51   const Char_t *AliMeanVertexPreprocessorOffline::fgkStatusCodeName[AliMeanVertexPreprocessorOffline::kNStatusCodes] = {
  52   "ok",
  53   "open file error or missing histos",
  54   "too low statistics",
  55   "problems storing OCDB",
  56   "write MeanVertex computed online",
  57   "write SPD vtx offline",
  58   "lumi region or cov matrix computation problems, default values set"
  59 };
  60
  61
  62 //____________________________________________________
  63 AliMeanVertexPreprocessorOffline::AliMeanVertexPreprocessorOffline():
  64   TNamed("AliMeanVertexPreprocessorOffline","AliMeanVertexPreprocessorOffline"),
  65   fStatus(kOk)
  66 {
  67   //constructor
  68 }
  69 //____________________________________________________
  70
  71 AliMeanVertexPreprocessorOffline::~AliMeanVertexPreprocessorOffline()
  72 {
  73   //destructor
  74
  75 }
  76 //____________________________________________________
  77 void AliMeanVertexPreprocessorOffline::ProcessOutput(const char *filename, AliCDBStorage *db, Int_t runNb){
  78
  79         TFile *file = TFile::Open(filename);
  80         if (!file || !file->IsOpen()){
  81                 AliError(Form("cannot open outputfile %s", filename));
  82                 fStatus=kInputError;
  83                 return;
  84         }
  85
  86         if (!db){
  87           AliError(Form("no OCDB storage found, return"));
  88           fStatus=kInputError;
  89           return;
  90         }
  91
  92
  93         TList *list = (TList*)file->Get("MeanVertex");
  94
  95         TH1F *histTRKvtxX = 0x0;
  96         TH1F *histTRKvtxY = 0x0;
  97         TH1F *histTRKvtxZ = 0x0;
  98
  99         TH1F *histSPDvtxX = 0x0;
 100         TH1F *histSPDvtxY = 0x0;
 101         TH1F *histSPDvtxZ = 0x0;
 102
 103
 104         Bool_t useTRKvtx = kTRUE;
 105         Bool_t useITSSAvtx = kFALSE;
 106         Bool_t useSPDvtx = kFALSE;
 107         Bool_t spdAvailable = kTRUE;
 108         Bool_t writeMeanVertexSPD = kFALSE;
 109         Bool_t vertexerSPD3Doff=kFALSE;
 110
 111
 112     if (!list) {
 113
 114       histTRKvtxX = (TH1F*)file->Get("hTRKVertexX");
 115       histTRKvtxY = (TH1F*)file->Get("hTRKVertexY");
 116       histTRKvtxZ = (TH1F*)file->Get("hTRKVertexZ");
 117
 118       histSPDvtxX = (TH1F*)file->Get("hSPDVertexX");
 119       histSPDvtxY = (TH1F*)file->Get("hSPDVertexY");
 120       histSPDvtxZ = (TH1F*)file->Get("hSPDVertexZ");
 121
 122       if (!histTRKvtxX || !histTRKvtxY || !histTRKvtxZ){
 123
 124         useTRKvtx = kFALSE;
 125         useITSSAvtx = kTRUE;
 126
 127         histTRKvtxX = (TH1F*)file->FindObject("hITSSAVertexX");
 128         histTRKvtxY = (TH1F*)file->FindObject("hITSSAVertexY");
 129         histTRKvtxZ = (TH1F*)file->FindObject("hITSSAVertexZ");
 130
 131         if (!histTRKvtxX || !histTRKvtxY || !histTRKvtxZ){
 132
 133           useITSSAvtx=kFALSE;
 134           useSPDvtx=kTRUE;
 135
 136           if (!histSPDvtxX || !histSPDvtxY || !histSPDvtxZ){
 137             AliError(Form("cannot find any histograms available from file"));
 138             fStatus=kInputError;
 139             return;
 140           }
 141         }
 142       }
 143     }
 144
 145     else{
 146
 147       histTRKvtxX = (TH1F*)list->FindObject("hTRKVertexX");
 148       histTRKvtxY = (TH1F*)list->FindObject("hTRKVertexY");
 149       histTRKvtxZ = (TH1F*)list->FindObject("hTRKVertexZ");
 150
 151       histSPDvtxX = (TH1F*)list->FindObject("hSPDVertexX");
 152       histSPDvtxY = (TH1F*)list->FindObject("hSPDVertexY");
 153       histSPDvtxZ = (TH1F*)list->FindObject("hSPDVertexZ");
 154
 155       if (!histTRKvtxX || !histTRKvtxY || !histTRKvtxZ){
 156
 157         useTRKvtx = kFALSE;
 158         useITSSAvtx = kTRUE;
 159
 160         histTRKvtxX = (TH1F*)list->FindObject("hITSSAVertexX");
 161         histTRKvtxY = (TH1F*)list->FindObject("hITSSAVertexY");
 162         histTRKvtxZ = (TH1F*)list->FindObject("hITSSAVertexZ");
 163
 164         if (!histTRKvtxX || !histTRKvtxY || !histTRKvtxZ){
 165
 166           useITSSAvtx=kFALSE;
 167           useSPDvtx=kTRUE;
 168
 169           if (!histSPDvtxX || !histSPDvtxY || !histSPDvtxZ){
 170             AliError(Form("cannot find any histograms available from list"));
 171             fStatus=kInputError;
 172             return;
 173           }
 174         }
 175
 176       }
 177     }
 178
 179
 180     if (useTRKvtx){
 181
 182       Float_t nEntriesX = histTRKvtxX->GetEffectiveEntries();
 183       Float_t nEntriesY = histTRKvtxY->GetEffectiveEntries();
 184       Float_t nEntriesZ = histTRKvtxZ->GetEffectiveEntries();
 185
 186       if (nEntriesX < 50. || nEntriesY<50. || nEntriesZ<50.) {
 187         AliError(Form("TRK vertex histograms have too few entries for fitting"));
 188         useTRKvtx=kFALSE;
 189         useSPDvtx = kTRUE;
 190       }
 191     }
 192     if (useITSSAvtx){
 193
 194       Float_t nEntriesX = histTRKvtxX->GetEffectiveEntries();
 195       Float_t nEntriesY = histTRKvtxY->GetEffectiveEntries();
 196       Float_t nEntriesZ = histTRKvtxZ->GetEffectiveEntries();
 197
 198       if (nEntriesX < 50. || nEntriesY<50. || nEntriesZ<50.) {
 199         AliError(Form("ITSSA vertex histograms have too few entries for fitting"));
 200         useITSSAvtx=kFALSE;
 201         useSPDvtx=kTRUE;
 202       }
 203     }
 204
 205     Float_t nEntriesX = histSPDvtxX->GetEffectiveEntries();
 206     Float_t nEntriesY = histSPDvtxY->GetEffectiveEntries();
 207     Float_t nEntriesZ = histSPDvtxZ->GetEffectiveEntries();
 208
 209     if (nEntriesX < 50. || nEntriesY<50. || nEntriesZ<50.) {
 210       spdAvailable = kFALSE;
 211       if ((useTRKvtx==kFALSE) && (useITSSAvtx==kFALSE)){
 212         AliError(Form("Also SPD vertex histograms have too few entries for fitting, return"));
 213         fStatus=kLowStatistics;
 214         return;
 215       }
 216     }
 217
 218     if((nEntriesX == 0.)&&(nEntriesY==0.) && (nEntriesZ>0.)){
 219       vertexerSPD3Doff = kTRUE;
 220       AliWarning("Vertexer SPD 3D off");
 221     }
 222
 223     Double_t xMeanVtx=0., yMeanVtx=0., zMeanVtx=0.;
 224     Double_t xSigmaVtx=0., ySigmaVtx=0., zSigmaVtx=0.;
 225
 226
 227     TF1 *fitVtxX, *fitVtxY, *fitVtxZ;
 228
 229     if (useTRKvtx || useITSSAvtx){
 230       histTRKvtxX ->Fit("gaus", "M");
 231       fitVtxX = histTRKvtxX -> GetFunction("gaus");
 232       xMeanVtx = fitVtxX -> GetParameter(1);
 233       if (TMath::Abs(xMeanVtx) > 2.) {
 234         xMeanVtx = 0.;
 235         writeMeanVertexSPD=kTRUE;
 236         fStatus=kWriteMeanVertexSPD;
 237       }
 238
 239       histTRKvtxY ->Fit("gaus", "M");
 240       fitVtxY = histTRKvtxY -> GetFunction("gaus");
 241       yMeanVtx = fitVtxY -> GetParameter(1);
 242       if (TMath::Abs(yMeanVtx) > 2.) {
 243         yMeanVtx = 0.;
 244         writeMeanVertexSPD=kTRUE;
 245         fStatus=kWriteMeanVertexSPD;
 246       }
 247
 248       histTRKvtxZ ->Fit("gaus", "M", "", -12, 12);
 249       fitVtxZ = histTRKvtxZ -> GetFunction("gaus");
 250       zMeanVtx = fitVtxZ -> GetParameter(1);
 251       zSigmaVtx = fitVtxZ -> GetParameter(2);
 252       if ((TMath::Abs(zMeanVtx) > 20.) || (zSigmaVtx>12.)) {
 253         zMeanVtx = histTRKvtxZ->GetMean();
 254         zSigmaVtx = histTRKvtxZ->GetRMS();
 255         writeMeanVertexSPD=kTRUE;
 256         fStatus=kWriteMeanVertexSPD;
 257       }
 258
 259     }
 260
 261
 262     //check fits: compare histo mean with fit mean value
 263     Double_t xHistoMean, yHistoMean, zHistoMean;
 264     Double_t xHistoRMS, yHistoRMS, zHistoRMS;
 265
 266     if (useTRKvtx || useITSSAvtx){
 267       xHistoMean = histTRKvtxX -> GetMean();
 268       xHistoRMS = histTRKvtxX ->GetRMS();
 269
 270       if ((TMath::Abs(xHistoMean-xMeanVtx) > 0.5)){
 271         AliWarning(Form("Possible problems with the fit mean very different from histo mean... using SPD vertex"));
 272         useTRKvtx = kFALSE;
 273         useITSSAvtx = kFALSE;
 274         useSPDvtx = kTRUE;
 275         fStatus=kUseOfflineSPDvtx;
 276       }
 277
 278       yHistoMean = histTRKvtxY ->GetMean();
 279       yHistoRMS = histTRKvtxY ->GetRMS();
 280
 281       if ((TMath::Abs(yHistoMean-yMeanVtx) > 0.5)){
 282         AliWarning(Form("Possible problems with the fit mean very different from histo mean... using SPD vertex"));
 283         useTRKvtx = kFALSE;
 284         useITSSAvtx = kFALSE;
 285         useSPDvtx = kTRUE;
 286         fStatus=kUseOfflineSPDvtx;
 287       }
 288
 289       zHistoMean = histTRKvtxZ -> GetMean();
 290       zHistoRMS = histTRKvtxZ ->GetRMS();
 291
 292       if ((TMath::Abs(zHistoMean-zMeanVtx) > 1.)){
 293         AliWarning(Form("Possible problems with the fit mean very different from histo mean... using SPD vertex"));
 294         useTRKvtx = kFALSE;
 295         useITSSAvtx = kFALSE;
 296         useSPDvtx = kTRUE;
 297         fStatus=kUseOfflineSPDvtx;
 298       }
 299       AliDebug(2, Form("xHistoRMS = %f, yHistoRMS = %f, zHistoRMS = %f", xHistoRMS, yHistoRMS, zHistoRMS));
 300     }
 301
 302
 303     if ((useSPDvtx) && (spdAvailable) && (!vertexerSPD3Doff)){
 304
 305       histSPDvtxX ->Fit("gaus", "M");
 306       fitVtxX = histSPDvtxX -> GetFunction("gaus");
 307       xMeanVtx = fitVtxX -> GetParameter(1);
 308       xSigmaVtx = fitVtxX -> GetParameter(2);
 309       if (TMath::Abs(xMeanVtx) > 2.) {
 310         xMeanVtx = 0.;
 311         writeMeanVertexSPD=kTRUE;
 312       }
 313
 314       histSPDvtxY ->Fit("gaus", "M");
 315       fitVtxY = histSPDvtxY -> GetFunction("gaus");
 316       yMeanVtx = fitVtxY -> GetParameter(1);
 317       ySigmaVtx = fitVtxY -> GetParameter(2);
 318       if (TMath::Abs(yMeanVtx) > 2.) {
 319         yMeanVtx = 0.;
 320         writeMeanVertexSPD=kTRUE;
 321       }
 322
 323       histSPDvtxZ ->Fit("gaus", "M", "", -12, 12);
 324       fitVtxZ = histSPDvtxZ -> GetFunction("gaus");
 325       zMeanVtx = fitVtxZ -> GetParameter(1);
 326       zSigmaVtx = fitVtxZ -> GetParameter(2);
 327       if ((TMath::Abs(zMeanVtx) > 20.) || (zSigmaVtx>12.)) {
 328         zMeanVtx = histSPDvtxZ ->GetMean();
 329         zSigmaVtx = histSPDvtxZ->GetRMS();
 330         writeMeanVertexSPD = kTRUE;
 331       }
 332
 333     }
 334     else if ((useSPDvtx) && (!spdAvailable)){
 335       AliError(Form("Difference between trkVtx and online one, SPD histos not enough entry or SPD 3D vertex off. Writing Mean Vertex SPD"));
 336       writeMeanVertexSPD = kTRUE;
 337     }
 338
 339
 340     //check with online position
 341
 342     Double_t posOnline[3], sigmaOnline[3];
 343
 344     if (useTRKvtx || useITSSAvtx || writeMeanVertexSPD){
 345       AliCDBManager *manCheck = AliCDBManager::Instance();
 346       manCheck->SetDefaultStorage("raw://");
 347       manCheck->SetRun(runNb);
 348
 349       AliCDBEntry *entr = manCheck->Get("GRP/Calib/MeanVertexSPD");
 350       if(entr) {
 351         AliESDVertex *vtxOnline = (AliESDVertex*)entr->GetObject();
 352
 353         posOnline[0] = vtxOnline->GetX();
 354         posOnline[1] = vtxOnline->GetY();
 355         posOnline[2] = vtxOnline->GetZ();
 356
 357         sigmaOnline[0] = vtxOnline->GetXRes();
 358         sigmaOnline[1] = vtxOnline->GetYRes();
 359         sigmaOnline[2] = vtxOnline->GetZRes();
 360
 361         AliDebug(2, Form("sigmaOnline[0] = %f, sigmaOnline[1] = %f, sigmaOnline[2] = %f", sigmaOnline[0], sigmaOnline[1], sigmaOnline[2]));
 362         //vtxOnline->GetSigmaXYZ(sigmaOnline);
 363
 364         if ((TMath::Abs(posOnline[0]-xMeanVtx) > 0.1) || (TMath::Abs(posOnline[1]-yMeanVtx) > 0.1) || (TMath::Abs(posOnline[2]-zMeanVtx) > 1.)){
 365           AliWarning(Form("vertex offline far from the online one"));
 366         }
 367       }
 368     }
 369
 370
 371
 372     if (writeMeanVertexSPD){
 373
 374       AliWarning(Form("Writing Mean Vertex SPD, Mean Vertex not available"));
 375
 376       Double_t sigma[3]={0.0150, 0.0150, zSigmaVtx};
 377
 378       AliESDVertex  *vertex =  new AliESDVertex(posOnline, sigma, "vertex");
 379
 380       AliCDBId id("GRP/Calib/MeanVertex", runNb, runNb);
 381
 382       AliCDBMetaData metaData;
 383       metaData.SetBeamPeriod(0); //check!!!!
 384       metaData.SetResponsible("Davide Caffarri");
 385       metaData.SetComment("Mean Vertex object used in reconstruction");
 386
 387       if (!db->Put(vertex, id, &metaData)) {
 388         AliError(Form("Error while putting object in storage %s", db->GetURI().Data()));
 389       }
 390
 391       delete vertex;
 392       return;
 393     }
 394
 395
 396     Float_t meanMult = 40.;
 397     Float_t p2 = 1.4;
 398     Float_t resolVtx = 0.05;
 399
 400     Double_t xSigmaMult, ySigmaMult, corrXZ, corrYZ, lumiRegSquaredX, lumiRegSquaredY;
 401     Double_t covarXZ=0., covarYZ=0.;
 402
 403     Bool_t highMultEnvironment = kFALSE;
 404     Bool_t highMultppEnvironment = kFALSE;
 405     Bool_t lowMultppEnvironment = kFALSE;
 406
 407     TF1 *sigmaFitX, *sigmaFitY, *corrFit;
 408
 409     //TH1F *histTRKdefMultX=0;
 410     //TH1F *histTRKdefMultY=0;
 411     TH1F *histTRKHighMultX=0;
 412     TH1F *histTRKHighMultY=0;
 413     TH2F *histTRKVertexXZ=0;
 414     TH2F *histTRKVertexYZ=0;
 415
 416     TH2F *histTRKvsMultX=0x0;
 417     TH2F *histTRKvsMultY=0x0;
 418
 419     if (useTRKvtx){
 420       if (list){
 421               //histTRKdefMultX = (TH1F*)list->FindObject("hTRKVertexXdefMult");
 422               //histTRKdefMultY = (TH1F*)list->FindObject("hTRKVertexYdefMult");
 423         histTRKHighMultX = (TH1F*)list->FindObject("hTRKVertexXHighMult");
 424         histTRKHighMultY = (TH1F*)list->FindObject("hTRKVertexYHighMult");
 425
 426         histTRKvsMultX = (TH2F*)list->FindObject("hTRKVertexXvsMult");
 427         histTRKvsMultY = (TH2F*)list->FindObject("hTRKVertexYvsMult");
 428
 429         histTRKVertexXZ = (TH2F*)list->FindObject("hTRKVertexXZ");
 430         histTRKVertexYZ = (TH2F*)list->FindObject("hTRKVertexYZ");
 431       }
 432
 433       else {
 434               //histTRKdefMultX = (TH1F*)file->Get("hTRKVertexXdefMult");
 435               //histTRKdefMultY = (TH1F*)file->Get("hTRKVertexYdefMult");
 436         histTRKHighMultX = (TH1F*)file->Get("hTRKVertexXHighMult");
 437         histTRKHighMultY = (TH1F*)file->Get("hTRKVertexYHighMult");
 438
 439         histTRKvsMultX = (TH2F*)file->FindObject("hTRKVertexXvsMult");
 440         histTRKvsMultY = (TH2F*)file->FindObject("hTRKVertexYvsMult");
 441
 442         histTRKVertexXZ = (TH2F*)file->Get("hTRKVertexXZ");
 443         histTRKVertexYZ = (TH2F*)file->Get("hTRKVertexYZ");
 444       }
 445
 446     }
 447
 448     if (useITSSAvtx){
 449       if (list){
 450               //histTRKdefMultX = (TH1F*)list->FindObject("hITSSAVertexXdefMult");
 451               //histTRKdefMultY = (TH1F*)list->FindObject("hITSSAVertexYdefMult");
 452         histTRKHighMultX = (TH1F*)list->FindObject("hITSSAVertexXHighMult");
 453         histTRKHighMultY = (TH1F*)list->FindObject("hITSSAVertexYHighMult");
 454
 455         histTRKvsMultX = (TH2F*)file->FindObject("hITSSAVertexXvsMult");
 456         histTRKvsMultY = (TH2F*)file->FindObject("hITSSAVertexYvsMult");
 457
 458         histTRKVertexXZ = (TH2F*)list->FindObject("hITSSAVertexXZ");
 459         histTRKVertexYZ = (TH2F*)list->FindObject("hITSSAVertexYZ");
 460       }
 461
 462       else {
 463               //histTRKdefMultX = (TH1F*)file->Get("hITSSAVertexXdefMult");
 464               //histTRKdefMultY = (TH1F*)file->Get("hITSSAVertexYdefMult");
 465         histTRKHighMultX = (TH1F*)file->Get("hITSSAVertexXHighMult");
 466         histTRKHighMultY = (TH1F*)file->Get("hITSSAVertexYHighMult");
 467
 468         histTRKvsMultX = (TH2F*)file->FindObject("hITSSAVertexXvsMult");
 469         histTRKvsMultY = (TH2F*)file->FindObject("hITSSAVertexYvsMult");
 470
 471         histTRKVertexXZ = (TH2F*)file->Get("hITSSAVertexXZ");
 472         histTRKVertexYZ = (TH2F*)file->Get("hITSSAVertexYZ");
 473       }
 474     }
 475
 476
 477     TH1D *projXvsMult;
 478     TH1D *projYvsMult;
 479
 480     Float_t nEntriesMultX=0, nEntriesMultY=0.;
 481
 482     if ((histTRKHighMultX) && (histTRKHighMultY)){
 483
 484       projXvsMult = (TH1D*)histTRKvsMultX->ProjectionX("projXHighMultPbPb", 150, 300);
 485       projYvsMult = (TH1D*)histTRKvsMultY->ProjectionX("projYHighMultPbPb", 150, 300);
 486
 487       nEntriesMultX = projXvsMult->GetEffectiveEntries();
 488       nEntriesMultY = projYvsMult->GetEffectiveEntries();
 489
 490       if ((nEntriesMultX >100) && (nEntriesMultY>100)) {
 491         AliWarning(Form("Setting High Mulitplicity environment"));
 492         highMultEnvironment = kTRUE;
 493       }
 494       else {
 495
 496         projXvsMult = (TH1D*)histTRKvsMultX->ProjectionX("projXHighMultPbPb", 10, 30);
 497         projYvsMult = (TH1D*)histTRKvsMultY->ProjectionX("projYHighMultPbPb", 10, 30);
 498
 499         nEntriesMultX = projXvsMult->GetEffectiveEntries();
 500         nEntriesMultY = projYvsMult->GetEffectiveEntries();
 501
 502         if ((nEntriesMultX >100) && (nEntriesMultY>100)) {
 503           AliWarning(Form("Setting high pp Mulitplicity environment or p-A high multiplicity"));
 504           highMultppEnvironment=kTRUE;
 505         }
 506         else {
 507
 508           projXvsMult = (TH1D*)histTRKvsMultX->ProjectionX("projXHighMultPbPb", 3, 5);
 509           projYvsMult = (TH1D*)histTRKvsMultY->ProjectionX("projYHighMultPbPb", 3, 5);
 510
 511           nEntriesMultX = projXvsMult->GetEffectiveEntries();
 512           nEntriesMultY = projYvsMult->GetEffectiveEntries();
 513
 514           if ((nEntriesMultX >100) && (nEntriesMultY>100)) {
 515             AliWarning(Form("Setting low pp Mulitplicity environment"));
 516             lowMultppEnvironment=kTRUE;
 517           }
 518
 519         }
 520
 521       }
 522
 523
 524       if (lowMultppEnvironment==kTRUE) {
 525
 526         if ((projXvsMult->GetEntries() < 40.) || (projYvsMult->GetEntries() < 40.)){
 527           AliWarning(Form("histos for lumi reg calculation not found, default value set"));
 528           xSigmaVtx=0.0120;
 529           ySigmaVtx=0.0120;
 530           fStatus=kLumiRegCovMatrixProblem;
 531         } else {
 532
 533           projXvsMult -> Fit("gaus", "M", "", -0.4, 0.4);
 534           sigmaFitX = projXvsMult -> GetFunction("gaus");
 535           xSigmaMult = sigmaFitX->GetParameter(2);
 536
 537           lumiRegSquaredX = (xSigmaMult*xSigmaMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
 538           if (lumiRegSquaredX < 0 || lumiRegSquaredX < 1E-5) {
 539             AliWarning(Form("Difficult luminous region determination X, keep convoluted sigma"));
 540             xSigmaVtx = xSigmaMult;
 541             fStatus=kLumiRegCovMatrixProblem;
 542           }
 543
 544           else {
 545             if (lumiRegSquaredX > 0 && lumiRegSquaredX < 0.0005){
 546               xSigmaVtx = TMath::Sqrt(lumiRegSquaredX);
 547               xSigmaVtx = xSigmaVtx*1.1;
 548             }
 549
 550             else{
 551               AliWarning(Form("Not possible to define a luminous region X. Default values set"));
 552               xSigmaVtx = 0.0120;
 553               fStatus=kLumiRegCovMatrixProblem;
 554             }
 555           }
 556
 557           projYvsMult -> Fit("gaus", "M", "", -0.2, 0.6);
 558           sigmaFitY = projYvsMult -> GetFunction("gaus");
 559           ySigmaMult = sigmaFitY->GetParameter(2);
 560
 561           lumiRegSquaredY= (ySigmaMult*ySigmaMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
 562
 563           if (lumiRegSquaredY < 0 || lumiRegSquaredY < 1E-5) {
 564             AliWarning(Form("Difficult luminous region determination Y, keep convoluted sigma"));
 565             ySigmaVtx = ySigmaMult;
 566             fStatus=kLumiRegCovMatrixProblem;
 567           }
 568
 569           else{
 570             if (lumiRegSquaredY > 0 && lumiRegSquaredY < 0.0005){
 571               ySigmaVtx = TMath::Sqrt(lumiRegSquaredY);
 572               ySigmaVtx = ySigmaVtx*1.1;
 573             }
 574             else{
 575               AliWarning(Form("Not possible to define a luminous region Y. Default values set"));
 576               ySigmaVtx = 0.0120;
 577               fStatus=kLumiRegCovMatrixProblem;
 578             }
 579           }
 580         }
 581
 582         TProfile *htrkXZ = histTRKVertexXZ ->ProfileY();
 583         htrkXZ -> Fit("pol1", "M", "", -10., 10.);
 584         corrFit = htrkXZ->GetFunction("pol1");
 585         corrXZ = corrFit->GetParameter(1);
 586
 587         if (TMath::Abs(corrXZ) > 0.01) {
 588           AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 589           corrXZ =0.;
 590           fStatus=kLumiRegCovMatrixProblem;
 591         }
 592         else{
 593           covarXZ = corrXZ * zSigmaVtx*zSigmaVtx;
 594
 595         }
 596
 597         TProfile *htrkYZ = histTRKVertexYZ ->ProfileY();
 598         htrkYZ -> Fit("pol1", "M", "", -10., 10.);
 599         corrFit = htrkYZ->GetFunction("pol1");
 600         corrYZ = corrFit->GetParameter(1);
 601
 602         if (TMath::Abs(corrYZ) > 0.01) {
 603           AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 604           corrYZ =0.;
 605           fStatus=kLumiRegCovMatrixProblem;
 606         }
 607         else{
 608           covarYZ = corrYZ*zSigmaVtx*zSigmaVtx;
 609         }
 610
 611       }
 612     }
 613
 614      if (highMultEnvironment==kTRUE || highMultppEnvironment==kTRUE){
 615
 616        projXvsMult -> Fit("gaus", "M", "", -0.4, 0.4);
 617        sigmaFitX = projXvsMult -> GetFunction("gaus");
 618        xSigmaMult = sigmaFitX->GetParameter(2);
 619
 620        if ((xSigmaMult <0) || (xSigmaMult>0.03)){
 621          AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
 622          xSigmaMult = 0.;
 623          xSigmaVtx = 0.0120;
 624          fStatus=kLumiRegCovMatrixProblem;
 625        }
 626        else{
 627          xSigmaVtx = xSigmaMult;
 628          xSigmaVtx = xSigmaVtx*1.1;
 629       }
 630
 631        projYvsMult -> Fit("gaus", "M", "", -0.2, 0.5);
 632        TCanvas *c = new TCanvas("nwC", "nwC");
 633        c->cd();
 634        projYvsMult->Draw();
 635       sigmaFitY = projYvsMult -> GetFunction("gaus");
 636       ySigmaMult = sigmaFitY->GetParameter(2);
 637
 638       if ((ySigmaMult <0) || (ySigmaMult>0.03)){
 639         AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
 640         ySigmaMult = 0.;
 641         ySigmaVtx = 0.0120;
 642         fStatus=kLumiRegCovMatrixProblem;
 643       }
 644       else{
 645         ySigmaVtx = ySigmaMult;
 646         ySigmaVtx = ySigmaVtx*1.1;
 647       }
 648
 649       TProfile *htrkXZ = histTRKVertexXZ ->ProfileY();
 650       htrkXZ -> Fit("pol1", "M", "", -10., 10.);
 651       corrFit = htrkXZ->GetFunction("pol1");
 652       corrXZ = corrFit->GetParameter(1);
 653
 654       if (TMath::Abs(corrXZ) > 0.01) {
 655         AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 656         corrXZ =0.;
 657         fStatus=kLumiRegCovMatrixProblem;
 658       }
 659       else{
 660         covarXZ = corrXZ * zSigmaVtx*zSigmaVtx;
 661
 662       }
 663
 664       TProfile *htrkYZ = histTRKVertexYZ ->ProfileY();
 665       htrkYZ -> Fit("pol1", "M", "", -10., 10.);
 666       corrFit = htrkYZ->GetFunction("pol1");
 667       corrYZ = corrFit->GetParameter(1);
 668
 669       if (TMath::Abs(corrYZ) > 0.01) {
 670         AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 671         corrYZ =0.;
 672         fStatus=kLumiRegCovMatrixProblem;
 673       }
 674       else{
 675         covarYZ = corrYZ*zSigmaVtx*zSigmaVtx;
 676       }
 677
 678      }
 679
 680      Double_t position[3], covMatrix[6];
 681      Double_t chi2=1.;
 682      Int_t nContr=1;
 683
 684      position[0] = xMeanVtx;
 685      position[1] = yMeanVtx;
 686      position[2] = zMeanVtx;
 687
 688      covMatrix[0] = xSigmaVtx*xSigmaVtx;
 689      covMatrix[1] = 0.; //xy
 690      covMatrix[2] = ySigmaVtx*ySigmaVtx;
 691      covMatrix[3] = covarXZ;
 692      covMatrix[4] = covarYZ;
 693      covMatrix[5] = zSigmaVtx*zSigmaVtx;
 694
 695
 696      //Printf ("%f, %f, %f, %f", xSigmaVtx, ySigmaVtx, covarXZ, covarYZ);
 697
 698      AliESDVertex  *vertex =  new AliESDVertex(position, covMatrix, chi2, nContr, "vertex");
 699
 700      AliCDBId id("GRP/Calib/MeanVertex", runNb, runNb);
 701
 702      AliCDBMetaData metaData;
 703      metaData.SetBeamPeriod(0); //check!!!!
 704      metaData.SetResponsible("Davide Caffarri");
 705      metaData.SetComment("Mean Vertex object used in reconstruction");
 706
 707      if (!db->Put(vertex, id, &metaData)) {
 708        AliError(Form("Error while putting object in storage %s", db->GetURI().Data()));
 709        fStatus=kStoreError;
 710      }
 711
 712      delete vertex;
 713
 714      Int_t status=GetStatus();
 715      if (status == 0) {
 716        AliInfo(Form("MeanVertex calibration successful: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 717      }
 718      else if (status > 0) {
 719        AliInfo(Form("MeanVertex calibration failed: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 720      }
 721      else if (status < 0) {
 722        AliInfo(Form("MeanVertex calibration but not fatal error: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 723      }
 724
 725
 726 }
 727
 728 //__________________________________________________________________________
 729 Int_t AliMeanVertexPreprocessorOffline::GetStatus(){
 730   /*
 731    * get status
 732    */
 733
 734   switch (fStatus) {
 735
 736     /* OK, return zero */
 737   case kOk:
 738     return 0;
 739     break;
 740
 741     /* non-fatal error, return negative status */
 742   case kLowStatistics:
 743   case kWriteMeanVertexSPD:
 744   case kUseOfflineSPDvtx:
 745   case kLumiRegCovMatrixProblem:
 746     return -fStatus;
 747     break;
 748
 749     /* fatal error, return positive status */
 750   case kInputError:
 751   case kStoreError:
 752     return fStatus;
 753     break;
 754
 755     /* anything else, return negative large number */
 756   default:
 757     return -999;
 758     break;
 759   }
 760
 761   /* should never arrive here, anyway return negative large number */
 762   return -999;
 763 }
 764