]> 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
added the alien caching script with example config
[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;
214 return;
215 }
216 }
f85d4cb9 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 }
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);
e0ea2e34 345
791a58dc 346 AliCDBEntry *entr = manCheck->Get("GRP/Calib/MeanVertexSPD");
347 if(entr) {
348 AliESDVertex *vtxOnline = (AliESDVertex*)entr->GetObject();
e0ea2e34 349
791a58dc 350 posOnline[0] = vtxOnline->GetX();
351 posOnline[1] = vtxOnline->GetY();
352 posOnline[2] = vtxOnline->GetZ();
e0ea2e34 353
791a58dc 354 sigmaOnline[0] = vtxOnline->GetXRes();
355 sigmaOnline[1] = vtxOnline->GetYRes();
356 sigmaOnline[2] = vtxOnline->GetZRes();
e0ea2e34 357
791a58dc 358 //vtxOnline->GetSigmaXYZ(sigmaOnline);
be5ba8df 359
791a58dc 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"));
e0ea2e34 362 }
791a58dc 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
ba8ff934 405 //TH1F *histTRKdefMultX=0;
406 //TH1F *histTRKdefMultY=0;
791a58dc 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){
f906054a 417 //histTRKdefMultX = (TH1F*)list->FindObject("hTRKVertexXdefMult");
418 //histTRKdefMultY = (TH1F*)list->FindObject("hTRKVertexYdefMult");
791a58dc 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 {
f906054a 430 //histTRKdefMultX = (TH1F*)file->Get("hTRKVertexXdefMult");
431 //histTRKdefMultY = (TH1F*)file->Get("hTRKVertexYdefMult");
791a58dc 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){
c308bb4a 446 //histTRKdefMultX = (TH1F*)list->FindObject("hITSSAVertexXdefMult");
447 //histTRKdefMultY = (TH1F*)list->FindObject("hITSSAVertexYdefMult");
791a58dc 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 {
c308bb4a 459 //histTRKdefMultX = (TH1F*)file->Get("hITSSAVertexXdefMult");
460 //histTRKdefMultY = (TH1F*)file->Get("hITSSAVertexYdefMult");
791a58dc 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;
be5ba8df 501 }
791a58dc 502 else {
adca66d2 503
791a58dc 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;
792b5897 513 }
a680415a 514
a680415a 515 }
e0ea2e34 516
791a58dc 517 }
518
519
520 if (lowMultppEnvironment==kTRUE) {
e0ea2e34 521
791a58dc 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 {
adca66d2 528
791a58dc 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;
4d14c5b5 537 fStatus=kLumiRegCovMatrixProblem;
791a58dc 538 }
539
540 else {
4d14c5b5 541 if (lumiRegSquaredX > 0 && lumiRegSquaredX < 0.0005){
542 xSigmaVtx = TMath::Sqrt(lumiRegSquaredX);
791a58dc 543 xSigmaVtx = xSigmaVtx*1.1;
4d14c5b5 544 }
545
791a58dc 546 else{
547 AliWarning(Form("Not possible to define a luminous region X. Default values set"));
548 xSigmaVtx = 0.0120;
4d14c5b5 549 fStatus=kLumiRegCovMatrixProblem;
550 }
4d14c5b5 551 }
adca66d2 552
791a58dc 553 projYvsMult -> Fit("gaus", "M", "", -0.2, 0.6);
554 sigmaFitY = projYvsMult -> GetFunction("gaus");
4d14c5b5 555 ySigmaMult = sigmaFitY->GetParameter(2);
adca66d2 556
791a58dc 557 lumiRegSquaredY= (ySigmaMult*ySigmaMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
adca66d2 558
791a58dc 559 if (lumiRegSquaredY < 0 || lumiRegSquaredY < 1E-5) {
560 AliWarning(Form("Difficult luminous region determination Y, keep convoluted sigma"));
561 ySigmaVtx = ySigmaMult;
4d14c5b5 562 fStatus=kLumiRegCovMatrixProblem;
adca66d2 563 }
4d14c5b5 564
791a58dc 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 }
adca66d2 575 }
e0ea2e34 576 }
a680415a 577
791a58dc 578 TProfile *htrkXZ = histTRKVertexXZ ->ProfileY();
579 htrkXZ -> Fit("pol1", "M", "", -10., 10.);
580 corrFit = htrkXZ->GetFunction("pol1");
581 corrXZ = corrFit->GetParameter(1);
4d14c5b5 582
791a58dc 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
4d14c5b5 591 }
592
791a58dc 593 TProfile *htrkYZ = histTRKVertexYZ ->ProfileY();
594 htrkYZ -> Fit("pol1", "M", "", -10., 10.);
595 corrFit = htrkYZ->GetFunction("pol1");
596 corrYZ = corrFit->GetParameter(1);
4d14c5b5 597
791a58dc 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;
4d14c5b5 602 }
791a58dc 603 else{
604 covarYZ = corrYZ*zSigmaVtx*zSigmaVtx;
4d14c5b5 605 }
606
791a58dc 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
e0ea2e34 722}
723
4d14c5b5 724//__________________________________________________________________________
725Int_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}
e0ea2e34 760