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     }
 300
 301
 302     if ((useSPDvtx) && (spdAvailable) && (!vertexerSPD3Doff)){
 303
 304       histSPDvtxX ->Fit("gaus", "M");
 305       fitVtxX = histSPDvtxX -> GetFunction("gaus");
 306       xMeanVtx = fitVtxX -> GetParameter(1);
 307       if (TMath::Abs(xMeanVtx) > 2.) {
 308         xMeanVtx = 0.;
 309         writeMeanVertexSPD=kTRUE;
 310       }
 311
 312       histSPDvtxY ->Fit("gaus", "M");
 313       fitVtxY = histSPDvtxY -> GetFunction("gaus");
 314       yMeanVtx = fitVtxY -> GetParameter(1);
 315       if (TMath::Abs(yMeanVtx) > 2.) {
 316         yMeanVtx = 0.;
 317         writeMeanVertexSPD=kTRUE;
 318       }
 319
 320       histSPDvtxZ ->Fit("gaus", "M", "", -12, 12);
 321       fitVtxZ = histSPDvtxZ -> GetFunction("gaus");
 322       zMeanVtx = fitVtxZ -> GetParameter(1);
 323       zSigmaVtx = fitVtxZ -> GetParameter(2);
 324       if ((TMath::Abs(zMeanVtx) > 20.) || (zSigmaVtx>12.)) {
 325         zMeanVtx = histSPDvtxZ ->GetMean();
 326         zSigmaVtx = histSPDvtxZ->GetRMS();
 327         writeMeanVertexSPD = kTRUE;
 328       }
 329
 330     }
 331     else if ((useSPDvtx) && (!spdAvailable)){
 332       AliError(Form("Difference between trkVtx and online one, SPD histos not enough entry or SPD 3D vertex off. Writing Mean Vertex SPD"));
 333       writeMeanVertexSPD = kTRUE;
 334     }
 335
 336
 337     //check with online position
 338
 339     Double_t posOnline[3], sigmaOnline[3];
 340
 341     if (useTRKvtx || useITSSAvtx || writeMeanVertexSPD){
 342       AliCDBManager *manCheck = AliCDBManager::Instance();
 343       manCheck->SetDefaultStorage("raw://");
 344       manCheck->SetRun(runNb);
 345
 346       AliCDBEntry *entr = manCheck->Get("GRP/Calib/MeanVertexSPD");
 347       if(entr) {
 348         AliESDVertex *vtxOnline = (AliESDVertex*)entr->GetObject();
 349
 350         posOnline[0] = vtxOnline->GetX();
 351         posOnline[1] = vtxOnline->GetY();
 352         posOnline[2] = vtxOnline->GetZ();
 353
 354         sigmaOnline[0] = vtxOnline->GetXRes();
 355         sigmaOnline[1] = vtxOnline->GetYRes();
 356         sigmaOnline[2] = vtxOnline->GetZRes();
 357
 358         //vtxOnline->GetSigmaXYZ(sigmaOnline);
 359
 360         if ((TMath::Abs(posOnline[0]-xMeanVtx) > 0.1) || (TMath::Abs(posOnline[1]-yMeanVtx) > 0.1) || (TMath::Abs(posOnline[2]-zMeanVtx) > 1.)){
 361           AliWarning(Form("vertex offline far from the online one"));
 362         }
 363       }
 364     }
 365
 366
 367
 368     if (writeMeanVertexSPD){
 369
 370       AliWarning(Form("Writing Mean Vertex SPD, Mean Vertex not available"));
 371
 372       Double_t sigma[3]={0.0150, 0.0150, zSigmaVtx};
 373
 374       AliESDVertex  *vertex =  new AliESDVertex(posOnline, sigma, "vertex");
 375
 376       AliCDBId id("GRP/Calib/MeanVertex", runNb, runNb);
 377
 378       AliCDBMetaData metaData;
 379       metaData.SetBeamPeriod(0); //check!!!!
 380       metaData.SetResponsible("Davide Caffarri");
 381       metaData.SetComment("Mean Vertex object used in reconstruction");
 382
 383       if (!db->Put(vertex, id, &metaData)) {
 384         AliError(Form("Error while putting object in storage %s", db->GetURI().Data()));
 385       }
 386
 387       delete vertex;
 388       return;
 389     }
 390
 391
 392     Float_t meanMult = 40.;
 393     Float_t p2 = 1.4;
 394     Float_t resolVtx = 0.05;
 395
 396     Double_t xSigmaMult, ySigmaMult, corrXZ, corrYZ, lumiRegSquaredX, lumiRegSquaredY;
 397     Double_t covarXZ=0., covarYZ=0.;
 398
 399     Bool_t highMultEnvironment = kFALSE;
 400     Bool_t highMultppEnvironment = kFALSE;
 401     Bool_t lowMultppEnvironment = kFALSE;
 402
 403     TF1 *sigmaFitX, *sigmaFitY, *corrFit;
 404
 405     TH1F *histTRKdefMultX=0;
 406     TH1F *histTRKdefMultY=0;
 407     TH1F *histTRKHighMultX=0;
 408     TH1F *histTRKHighMultY=0;
 409     TH2F *histTRKVertexXZ=0;
 410     TH2F *histTRKVertexYZ=0;
 411
 412     TH2F *histTRKvsMultX=0x0;
 413     TH2F *histTRKvsMultY=0x0;
 414
 415     if (useTRKvtx){
 416       if (list){
 417               //histTRKdefMultX = (TH1F*)list->FindObject("hTRKVertexXdefMult");
 418               //histTRKdefMultY = (TH1F*)list->FindObject("hTRKVertexYdefMult");
 419         histTRKHighMultX = (TH1F*)list->FindObject("hTRKVertexXHighMult");
 420         histTRKHighMultY = (TH1F*)list->FindObject("hTRKVertexYHighMult");
 421
 422         histTRKvsMultX = (TH2F*)list->FindObject("hTRKVertexXvsMult");
 423         histTRKvsMultY = (TH2F*)list->FindObject("hTRKVertexYvsMult");
 424
 425         histTRKVertexXZ = (TH2F*)list->FindObject("hTRKVertexXZ");
 426         histTRKVertexYZ = (TH2F*)list->FindObject("hTRKVertexYZ");
 427       }
 428
 429       else {
 430               //histTRKdefMultX = (TH1F*)file->Get("hTRKVertexXdefMult");
 431               //histTRKdefMultY = (TH1F*)file->Get("hTRKVertexYdefMult");
 432         histTRKHighMultX = (TH1F*)file->Get("hTRKVertexXHighMult");
 433         histTRKHighMultY = (TH1F*)file->Get("hTRKVertexYHighMult");
 434
 435         histTRKvsMultX = (TH2F*)file->FindObject("hTRKVertexXvsMult");
 436         histTRKvsMultY = (TH2F*)file->FindObject("hTRKVertexYvsMult");
 437
 438         histTRKVertexXZ = (TH2F*)file->Get("hTRKVertexXZ");
 439         histTRKVertexYZ = (TH2F*)file->Get("hTRKVertexYZ");
 440       }
 441
 442     }
 443
 444     if (useITSSAvtx){
 445       if (list){
 446               //histTRKdefMultX = (TH1F*)list->FindObject("hITSSAVertexXdefMult");
 447               //histTRKdefMultY = (TH1F*)list->FindObject("hITSSAVertexYdefMult");
 448         histTRKHighMultX = (TH1F*)list->FindObject("hITSSAVertexXHighMult");
 449         histTRKHighMultY = (TH1F*)list->FindObject("hITSSAVertexYHighMult");
 450
 451         histTRKvsMultX = (TH2F*)file->FindObject("hITSSAVertexXvsMult");
 452         histTRKvsMultY = (TH2F*)file->FindObject("hITSSAVertexYvsMult");
 453
 454         histTRKVertexXZ = (TH2F*)list->FindObject("hITSSAVertexXZ");
 455         histTRKVertexYZ = (TH2F*)list->FindObject("hITSSAVertexYZ");
 456       }
 457
 458       else {
 459               //histTRKdefMultX = (TH1F*)file->Get("hITSSAVertexXdefMult");
 460               //histTRKdefMultY = (TH1F*)file->Get("hITSSAVertexYdefMult");
 461         histTRKHighMultX = (TH1F*)file->Get("hITSSAVertexXHighMult");
 462         histTRKHighMultY = (TH1F*)file->Get("hITSSAVertexYHighMult");
 463
 464         histTRKvsMultX = (TH2F*)file->FindObject("hITSSAVertexXvsMult");
 465         histTRKvsMultY = (TH2F*)file->FindObject("hITSSAVertexYvsMult");
 466
 467         histTRKVertexXZ = (TH2F*)file->Get("hITSSAVertexXZ");
 468         histTRKVertexYZ = (TH2F*)file->Get("hITSSAVertexYZ");
 469       }
 470     }
 471
 472
 473     TH1D *projXvsMult;
 474     TH1D *projYvsMult;
 475
 476     Float_t nEntriesMultX=0, nEntriesMultY=0.;
 477
 478     if ((histTRKHighMultX) && (histTRKHighMultY)){
 479
 480       projXvsMult = (TH1D*)histTRKvsMultX->ProjectionX("projXHighMultPbPb", 150, 300);
 481       projYvsMult = (TH1D*)histTRKvsMultY->ProjectionX("projYHighMultPbPb", 150, 300);
 482
 483       nEntriesMultX = projXvsMult->GetEffectiveEntries();
 484       nEntriesMultY = projYvsMult->GetEffectiveEntries();
 485
 486       if ((nEntriesMultX >100) && (nEntriesMultY>100)) {
 487         AliWarning(Form("Setting High Mulitplicity environment"));
 488         highMultEnvironment = kTRUE;
 489       }
 490       else {
 491
 492         projXvsMult = (TH1D*)histTRKvsMultX->ProjectionX("projXHighMultPbPb", 10, 30);
 493         projYvsMult = (TH1D*)histTRKvsMultY->ProjectionX("projYHighMultPbPb", 10, 30);
 494
 495         nEntriesMultX = projXvsMult->GetEffectiveEntries();
 496         nEntriesMultY = projYvsMult->GetEffectiveEntries();
 497
 498         if ((nEntriesMultX >100) && (nEntriesMultY>100)) {
 499           AliWarning(Form("Setting high pp Mulitplicity environment or p-A high multiplicity"));
 500           highMultppEnvironment=kTRUE;
 501         }
 502         else {
 503
 504           projXvsMult = (TH1D*)histTRKvsMultX->ProjectionX("projXHighMultPbPb", 3, 5);
 505           projYvsMult = (TH1D*)histTRKvsMultY->ProjectionX("projYHighMultPbPb", 3, 5);
 506
 507           nEntriesMultX = projXvsMult->GetEffectiveEntries();
 508           nEntriesMultY = projYvsMult->GetEffectiveEntries();
 509
 510           if ((nEntriesMultX >100) && (nEntriesMultY>100)) {
 511             AliWarning(Form("Setting low pp Mulitplicity environment"));
 512             lowMultppEnvironment=kTRUE;
 513           }
 514
 515         }
 516
 517       }
 518
 519
 520       if (lowMultppEnvironment==kTRUE) {
 521
 522         if ((projXvsMult->GetEntries() < 40.) || (projYvsMult->GetEntries() < 40.)){
 523           AliWarning(Form("histos for lumi reg calculation not found, default value set"));
 524           xSigmaVtx=0.0120;
 525           ySigmaVtx=0.0120;
 526           fStatus=kLumiRegCovMatrixProblem;
 527         } else {
 528
 529           projXvsMult -> Fit("gaus", "M", "", -0.4, 0.4);
 530           sigmaFitX = projXvsMult -> GetFunction("gaus");
 531           xSigmaMult = sigmaFitX->GetParameter(2);
 532
 533           lumiRegSquaredX = (xSigmaMult*xSigmaMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
 534           if (lumiRegSquaredX < 0 || lumiRegSquaredX < 1E-5) {
 535             AliWarning(Form("Difficult luminous region determination X, keep convoluted sigma"));
 536             xSigmaVtx = xSigmaMult;
 537             fStatus=kLumiRegCovMatrixProblem;
 538           }
 539
 540           else {
 541             if (lumiRegSquaredX > 0 && lumiRegSquaredX < 0.0005){
 542               xSigmaVtx = TMath::Sqrt(lumiRegSquaredX);
 543               xSigmaVtx = xSigmaVtx*1.1;
 544             }
 545
 546             else{
 547               AliWarning(Form("Not possible to define a luminous region X. Default values set"));
 548               xSigmaVtx = 0.0120;
 549               fStatus=kLumiRegCovMatrixProblem;
 550             }
 551           }
 552
 553           projYvsMult -> Fit("gaus", "M", "", -0.2, 0.6);
 554           sigmaFitY = projYvsMult -> GetFunction("gaus");
 555           ySigmaMult = sigmaFitY->GetParameter(2);
 556
 557           lumiRegSquaredY= (ySigmaMult*ySigmaMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
 558
 559           if (lumiRegSquaredY < 0 || lumiRegSquaredY < 1E-5) {
 560             AliWarning(Form("Difficult luminous region determination Y, keep convoluted sigma"));
 561             ySigmaVtx = ySigmaMult;
 562             fStatus=kLumiRegCovMatrixProblem;
 563           }
 564
 565           else{
 566             if (lumiRegSquaredY > 0 && lumiRegSquaredY < 0.0005){
 567               ySigmaVtx = TMath::Sqrt(lumiRegSquaredY);
 568               ySigmaVtx = ySigmaVtx*1.1;
 569             }
 570             else{
 571               AliWarning(Form("Not possible to define a luminous region Y. Default values set"));
 572               ySigmaVtx = 0.0120;
 573               fStatus=kLumiRegCovMatrixProblem;
 574             }
 575           }
 576         }
 577
 578         TProfile *htrkXZ = histTRKVertexXZ ->ProfileY();
 579         htrkXZ -> Fit("pol1", "M", "", -10., 10.);
 580         corrFit = htrkXZ->GetFunction("pol1");
 581         corrXZ = corrFit->GetParameter(1);
 582
 583         if (TMath::Abs(corrXZ) > 0.01) {
 584           AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 585           corrXZ =0.;
 586           fStatus=kLumiRegCovMatrixProblem;
 587         }
 588         else{
 589           covarXZ = corrXZ * zSigmaVtx*zSigmaVtx;
 590
 591         }
 592
 593         TProfile *htrkYZ = histTRKVertexYZ ->ProfileY();
 594         htrkYZ -> Fit("pol1", "M", "", -10., 10.);
 595         corrFit = htrkYZ->GetFunction("pol1");
 596         corrYZ = corrFit->GetParameter(1);
 597
 598         if (TMath::Abs(corrYZ) > 0.01) {
 599           AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 600           corrYZ =0.;
 601           fStatus=kLumiRegCovMatrixProblem;
 602         }
 603         else{
 604           covarYZ = corrYZ*zSigmaVtx*zSigmaVtx;
 605         }
 606
 607       }
 608     }
 609
 610      if (highMultEnvironment==kTRUE || highMultppEnvironment==kTRUE){
 611
 612     projXvsMult -> Fit("gaus", "M", "", -0.4, 0.4);
 613     sigmaFitX = projXvsMult -> GetFunction("gaus");
 614     xSigmaMult = sigmaFitX->GetParameter(2);
 615
 616       if ((xSigmaMult <0) || (xSigmaMult>0.03)){
 617         AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
 618         xSigmaMult = 0.;
 619         xSigmaVtx = 0.0120;
 620         fStatus=kLumiRegCovMatrixProblem;
 621       }
 622       else{
 623         xSigmaVtx = xSigmaMult;
 624         xSigmaVtx = xSigmaVtx*1.1;
 625       }
 626
 627       projYvsMult -> Fit("gaus", "M", "", -0.2, 0.5);
 628       TCanvas *c = new TCanvas("nwC", "nwC");
 629       c->cd();
 630       projYvsMult->Draw();
 631       sigmaFitY = projYvsMult -> GetFunction("gaus");
 632       ySigmaMult = sigmaFitY->GetParameter(2);
 633
 634       if ((ySigmaMult <0) || (ySigmaMult>0.03)){
 635         AliWarning(Form("Problems with luminosiy region determination, update of the postion only"));
 636         ySigmaMult = 0.;
 637         ySigmaVtx = 0.0120;
 638         fStatus=kLumiRegCovMatrixProblem;
 639       }
 640       else{
 641         ySigmaVtx = ySigmaMult;
 642         ySigmaVtx = ySigmaVtx*1.1;
 643       }
 644
 645       TProfile *htrkXZ = histTRKVertexXZ ->ProfileY();
 646       htrkXZ -> Fit("pol1", "M", "", -10., 10.);
 647       corrFit = htrkXZ->GetFunction("pol1");
 648       corrXZ = corrFit->GetParameter(1);
 649
 650       if (TMath::Abs(corrXZ) > 0.01) {
 651         AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 652         corrXZ =0.;
 653         fStatus=kLumiRegCovMatrixProblem;
 654       }
 655       else{
 656         covarXZ = corrXZ * zSigmaVtx*zSigmaVtx;
 657
 658       }
 659
 660       TProfile *htrkYZ = histTRKVertexYZ ->ProfileY();
 661       htrkYZ -> Fit("pol1", "M", "", -10., 10.);
 662       corrFit = htrkYZ->GetFunction("pol1");
 663       corrYZ = corrFit->GetParameter(1);
 664
 665       if (TMath::Abs(corrYZ) > 0.01) {
 666         AliWarning(Form("Problems in the correlation fitting, not update the covariance matrix"));
 667         corrYZ =0.;
 668         fStatus=kLumiRegCovMatrixProblem;
 669       }
 670       else{
 671         covarYZ = corrYZ*zSigmaVtx*zSigmaVtx;
 672       }
 673
 674      }
 675
 676      Double_t position[3], covMatrix[6];
 677      Double_t chi2=1.;
 678      Int_t nContr=1;
 679
 680      position[0] = xMeanVtx;
 681      position[1] = yMeanVtx;
 682      position[2] = zMeanVtx;
 683
 684      covMatrix[0] = xSigmaVtx*xSigmaVtx;
 685      covMatrix[1] = 0.; //xy
 686      covMatrix[2] = ySigmaVtx*ySigmaVtx;
 687      covMatrix[3] = covarXZ;
 688      covMatrix[4] = covarYZ;
 689      covMatrix[5] = zSigmaVtx*zSigmaVtx;
 690
 691
 692      //Printf ("%f, %f, %f, %f", xSigmaVtx, ySigmaVtx, covarXZ, covarYZ);
 693
 694      AliESDVertex  *vertex =  new AliESDVertex(position, covMatrix, chi2, nContr, "vertex");
 695
 696      AliCDBId id("GRP/Calib/MeanVertex", runNb, runNb);
 697
 698      AliCDBMetaData metaData;
 699      metaData.SetBeamPeriod(0); //check!!!!
 700      metaData.SetResponsible("Davide Caffarri");
 701      metaData.SetComment("Mean Vertex object used in reconstruction");
 702
 703      if (!db->Put(vertex, id, &metaData)) {
 704        AliError(Form("Error while putting object in storage %s", db->GetURI().Data()));
 705        fStatus=kStoreError;
 706      }
 707
 708      delete vertex;
 709
 710      Int_t status=GetStatus();
 711      if (status == 0) {
 712        AliInfo(Form("MeanVertex calibration successful: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 713      }
 714      else if (status > 0) {
 715        AliInfo(Form("MeanVertex calibration failed: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 716      }
 717      else if (status < 0) {
 718        AliInfo(Form("MeanVertex calibration but not fatal error: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 719      }
 720
 721
 722 }
 723
 724 //__________________________________________________________________________
 725 Int_t AliMeanVertexPreprocessorOffline::GetStatus(){
 726   /*
 727    * get status
 728    */
 729
 730   switch (fStatus) {
 731
 732     /* OK, return zero */
 733   case kOk:
 734     return 0;
 735     break;
 736
 737     /* non-fatal error, return negative status */
 738   case kLowStatistics:
 739   case kWriteMeanVertexSPD:
 740   case kUseOfflineSPDvtx:
 741   case kLumiRegCovMatrixProblem:
 742     return -fStatus;
 743     break;
 744
 745     /* fatal error, return positive status */
 746   case kInputError:
 747   case kStoreError:
 748     return fStatus;
 749     break;
 750
 751     /* anything else, return negative large number */
 752   default:
 753     return -999;
 754     break;
 755   }
 756
 757   /* should never arrive here, anyway return negative large number */
 758   return -999;
 759 }
 760