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