285a097db37172683e736f736d8b5c19d5a36f41
[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       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