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