]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/ITS/AliMeanVertexPreprocessorOffline.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGPP / ITS / AliMeanVertexPreprocessorOffline.cxx
CommitLineData
e0ea2e34 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"
791a58dc 38#include <TCanvas.h>
e0ea2e34 39
40#include "AliESDVertex.h"
41#include "TH1F.h"
42#include "TH2F.h"
43#include "TF1.h"
44#include "TProfile.h"
45
46
47ClassImp(AliMeanVertexPreprocessorOffline)
48
4d14c5b5 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
e0ea2e34 62//____________________________________________________
63AliMeanVertexPreprocessorOffline::AliMeanVertexPreprocessorOffline():
4d14c5b5 64 TNamed("AliMeanVertexPreprocessorOffline","AliMeanVertexPreprocessorOffline"),
65 fStatus(kOk)
e0ea2e34 66{
67 //constructor
68}
69//____________________________________________________
70
71AliMeanVertexPreprocessorOffline::~AliMeanVertexPreprocessorOffline()
72{
73 //destructor
74
75}
76//____________________________________________________
a9f9c69b 77void AliMeanVertexPreprocessorOffline::ProcessOutput(const char *filename, AliCDBStorage *db, Int_t runNb){
e0ea2e34 78
79 TFile *file = TFile::Open(filename);
80 if (!file || !file->IsOpen()){
81 AliError(Form("cannot open outputfile %s", filename));
4d14c5b5 82 fStatus=kInputError;
e0ea2e34 83 return;
84 }
85
a9f9c69b 86 if (!db){
87 AliError(Form("no OCDB storage found, return"));
4d14c5b5 88 fStatus=kInputError;
a9f9c69b 89 return;
e0ea2e34 90 }
792b5897 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
e0ea2e34 103
104 Bool_t useTRKvtx = kTRUE;
105 Bool_t useITSSAvtx = kFALSE;
106 Bool_t useSPDvtx = kFALSE;
be5ba8df 107 Bool_t spdAvailable = kTRUE;
108 Bool_t writeMeanVertexSPD = kFALSE;
a680415a 109 Bool_t vertexerSPD3Doff=kFALSE;
e0ea2e34 110
792b5897 111
112 if (!list) {
e0ea2e34 113
791a58dc 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){
adca66d2 132
791a58dc 133 useITSSAvtx=kFALSE;
134 useSPDvtx=kTRUE;
adca66d2 135
791a58dc 136 if (!histSPDvtxX || !histSPDvtxY || !histSPDvtxZ){
137 AliError(Form("cannot find any histograms available from file"));
138 fStatus=kInputError;
139 return;
140 }
adca66d2 141 }
791a58dc 142 }
143 }
144
145 else{
be5ba8df 146
791a58dc 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){
bbd90c24 165
791a58dc 166 useITSSAvtx=kFALSE;
167 useSPDvtx=kTRUE;
bbd90c24 168
791a58dc 169 if (!histSPDvtxX || !histSPDvtxY || !histSPDvtxZ){
170 AliError(Form("cannot find any histograms available from list"));
171 fStatus=kInputError;
172 return;
bbd90c24 173 }
791a58dc 174 }
e0ea2e34 175
791a58dc 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;
0e1b7412 214 return;
791a58dc 215 }
216 }
0e1b7412 217
791a58dc 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 }
fdc6afb7 299 AliDebug(2, Form("xHistoRMS = %f, yHistoRMS = %f, zHistoRMS = %f", xHistoRMS, yHistoRMS, zHistoRMS));
791a58dc 300 }
301
302
303 if ((useSPDvtx) && (spdAvailable) && (!vertexerSPD3Doff)){
304
305 histSPDvtxX ->Fit("gaus", "M");
306 fitVtxX = histSPDvtxX -> GetFunction("gaus");
307 xMeanVtx = fitVtxX -> GetParameter(1);
0e1b7412 308 xSigmaVtx = fitVtxX -> GetParameter(2);
791a58dc 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);
0e1b7412 317 ySigmaVtx = fitVtxY -> GetParameter(2);
791a58dc 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);
e0ea2e34 348
791a58dc 349 AliCDBEntry *entr = manCheck->Get("GRP/Calib/MeanVertexSPD");
350 if(entr) {
351 AliESDVertex *vtxOnline = (AliESDVertex*)entr->GetObject();
e0ea2e34 352
791a58dc 353 posOnline[0] = vtxOnline->GetX();
354 posOnline[1] = vtxOnline->GetY();
355 posOnline[2] = vtxOnline->GetZ();
e0ea2e34 356
791a58dc 357 sigmaOnline[0] = vtxOnline->GetXRes();
358 sigmaOnline[1] = vtxOnline->GetYRes();
359 sigmaOnline[2] = vtxOnline->GetZRes();
e0ea2e34 360
fdc6afb7 361 AliDebug(2, Form("sigmaOnline[0] = %f, sigmaOnline[1] = %f, sigmaOnline[2] = %f", sigmaOnline[0], sigmaOnline[1], sigmaOnline[2]));
791a58dc 362 //vtxOnline->GetSigmaXYZ(sigmaOnline);
be5ba8df 363
791a58dc 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"));
e0ea2e34 366 }
791a58dc 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
ba8ff934 409 //TH1F *histTRKdefMultX=0;
410 //TH1F *histTRKdefMultY=0;
791a58dc 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){
f906054a 421 //histTRKdefMultX = (TH1F*)list->FindObject("hTRKVertexXdefMult");
422 //histTRKdefMultY = (TH1F*)list->FindObject("hTRKVertexYdefMult");
791a58dc 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 {
f906054a 434 //histTRKdefMultX = (TH1F*)file->Get("hTRKVertexXdefMult");
435 //histTRKdefMultY = (TH1F*)file->Get("hTRKVertexYdefMult");
791a58dc 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){
c308bb4a 450 //histTRKdefMultX = (TH1F*)list->FindObject("hITSSAVertexXdefMult");
451 //histTRKdefMultY = (TH1F*)list->FindObject("hITSSAVertexYdefMult");
791a58dc 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 {
c308bb4a 463 //histTRKdefMultX = (TH1F*)file->Get("hITSSAVertexXdefMult");
464 //histTRKdefMultY = (TH1F*)file->Get("hITSSAVertexYdefMult");
791a58dc 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;
be5ba8df 505 }
791a58dc 506 else {
adca66d2 507
791a58dc 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;
792b5897 517 }
a680415a 518
a680415a 519 }
e0ea2e34 520
791a58dc 521 }
522
523
524 if (lowMultppEnvironment==kTRUE) {
e0ea2e34 525
791a58dc 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 {
adca66d2 532
791a58dc 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;
4d14c5b5 541 fStatus=kLumiRegCovMatrixProblem;
791a58dc 542 }
543
544 else {
4d14c5b5 545 if (lumiRegSquaredX > 0 && lumiRegSquaredX < 0.0005){
546 xSigmaVtx = TMath::Sqrt(lumiRegSquaredX);
791a58dc 547 xSigmaVtx = xSigmaVtx*1.1;
4d14c5b5 548 }
549
791a58dc 550 else{
551 AliWarning(Form("Not possible to define a luminous region X. Default values set"));
552 xSigmaVtx = 0.0120;
4d14c5b5 553 fStatus=kLumiRegCovMatrixProblem;
554 }
4d14c5b5 555 }
adca66d2 556
791a58dc 557 projYvsMult -> Fit("gaus", "M", "", -0.2, 0.6);
558 sigmaFitY = projYvsMult -> GetFunction("gaus");
4d14c5b5 559 ySigmaMult = sigmaFitY->GetParameter(2);
adca66d2 560
791a58dc 561 lumiRegSquaredY= (ySigmaMult*ySigmaMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
adca66d2 562
791a58dc 563 if (lumiRegSquaredY < 0 || lumiRegSquaredY < 1E-5) {
564 AliWarning(Form("Difficult luminous region determination Y, keep convoluted sigma"));
565 ySigmaVtx = ySigmaMult;
4d14c5b5 566 fStatus=kLumiRegCovMatrixProblem;
adca66d2 567 }
4d14c5b5 568
791a58dc 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 }
adca66d2 579 }
e0ea2e34 580 }
a680415a 581
791a58dc 582 TProfile *htrkXZ = histTRKVertexXZ ->ProfileY();
583 htrkXZ -> Fit("pol1", "M", "", -10., 10.);
584 corrFit = htrkXZ->GetFunction("pol1");
585 corrXZ = corrFit->GetParameter(1);
4d14c5b5 586
791a58dc 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
4d14c5b5 595 }
596
791a58dc 597 TProfile *htrkYZ = histTRKVertexYZ ->ProfileY();
598 htrkYZ -> Fit("pol1", "M", "", -10., 10.);
599 corrFit = htrkYZ->GetFunction("pol1");
600 corrYZ = corrFit->GetParameter(1);
4d14c5b5 601
791a58dc 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;
4d14c5b5 606 }
791a58dc 607 else{
608 covarYZ = corrYZ*zSigmaVtx*zSigmaVtx;
4d14c5b5 609 }
610
791a58dc 611 }
612 }
613
614 if (highMultEnvironment==kTRUE || highMultppEnvironment==kTRUE){
615
0e1b7412 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;
791a58dc 629 }
0e1b7412 630
631 projYvsMult -> Fit("gaus", "M", "", -0.2, 0.5);
632 TCanvas *c = new TCanvas("nwC", "nwC");
633 c->cd();
634 projYvsMult->Draw();
791a58dc 635 sigmaFitY = projYvsMult -> GetFunction("gaus");
636 ySigmaMult = sigmaFitY->GetParameter(2);
0e1b7412 637
791a58dc 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 }
0e1b7412 648
791a58dc 649 TProfile *htrkXZ = histTRKVertexXZ ->ProfileY();
650 htrkXZ -> Fit("pol1", "M", "", -10., 10.);
651 corrFit = htrkXZ->GetFunction("pol1");
652 corrXZ = corrFit->GetParameter(1);
0e1b7412 653
791a58dc 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 }
0e1b7412 663
791a58dc 664 TProfile *htrkYZ = histTRKVertexYZ ->ProfileY();
665 htrkYZ -> Fit("pol1", "M", "", -10., 10.);
666 corrFit = htrkYZ->GetFunction("pol1");
667 corrYZ = corrFit->GetParameter(1);
0e1b7412 668
791a58dc 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 }
0e1b7412 677
791a58dc 678 }
0e1b7412 679
791a58dc 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
e0ea2e34 726}
727
4d14c5b5 728//__________________________________________________________________________
729Int_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}
e0ea2e34 764