fine tuning of TOF tail (developing task)
[u/mrichter/AliRoot.git] / PWGPP / ITS / AliMeanVertexCalibTask.cxx
CommitLineData
bbb57828 1/**************************************************************************
2 * Copyright(c) 1998-2009, 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// Class AliMeanVertexCalibTask
18// AliAnalysisTask to extract from ESD the information on primary vertex
19// reconstruction in order to compute the MeanVertex object
20//
21// Author: D.Caffarri, davide.caffarri@pd.infn.it
22// A.Dainese, andrea.dainese@pd.infn.it
23//*************************************************************************
24
25
26#include <TH1F.h>
27#include <TH2F.h>
29be8a66 28#include <string.h>
bbb57828 29
30#include "AliMultiplicity.h"
31#include "AliESDVertex.h"
32#include "AliESDEvent.h"
33#include "AliVertexerTracks.h"
29be8a66 34#include "AliCDBManager.h"
35#include "AliCDBEntry.h"
36#include "AliGRPRecoParam.h"
bbb57828 37
38#include "AliMeanVertexCalibTask.h"
39
40
41ClassImp(AliMeanVertexCalibTask)
42
43//_____________________________________________________________________
44AliMeanVertexCalibTask::AliMeanVertexCalibTask(const char *name):
45 AliAnalysisTaskSE(name),
46 fESD(0),
47 fOutput(0),
48 fOnlyITSTPCTracks(kFALSE),
49 fOnlyITSSATracks(kTRUE)
50{
51
52 // Constructor
53
54 // Define input and output slots here
55 // Output slot #0 writes into a TList container
56 DefineOutput(1, TList::Class()); //My private output
57}
58
59//________________________________________________________________________
60AliMeanVertexCalibTask::~AliMeanVertexCalibTask()
61{
62 // Destructor
63
64 if (fOutput) {
65 delete fOutput;
66 fOutput = 0;
67 }
68
69}
70//________________________________________________________________________
71void AliMeanVertexCalibTask::UserCreateOutputObjects()
72{
73
74 // Several histograms are more conveniently managed in a TList
75 fOutput = new TList;
76 fOutput->SetOwner();
77
78
79 TH1F* hSPDVertexX = new TH1F("hSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1);
80 fOutput->Add(hSPDVertexX);
81 TH1F* hSPDVertexY = new TH1F("hSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1);
82 fOutput->Add(hSPDVertexY);
83 TH1F* hSPDVertexZ = new TH1F("hSPDVertexZ","SPDVertex z; z vertex [cm]; events",200,-20,20);
84 fOutput->Add(hSPDVertexZ);
85 TH1F* hTRKVertexX = new TH1F("hTRKVertexX","TRKVertex x; x vertex [cm]; events",200,-1,1);
86 fOutput->Add(hTRKVertexX);
87 TH1F* hTRKVertexY = new TH1F("hTRKVertexY","TRKVertex y; y vertex [cm]; events",200,-1,1);
88 fOutput->Add(hTRKVertexY);
89 TH1F* hTRKVertexZ = new TH1F("hTRKVertexZ","TRKVertex z; z vertex [cm]; events",200,-20,20);
90 fOutput->Add(hTRKVertexZ);
791a58dc 91
92 TH2F *hTRKVertexXvsMult = new TH2F("hTRKVertexXvsMult", "TRKVertex X vs mult", 200, -1, 1, 300, 0, 3000);
93 fOutput->Add(hTRKVertexXvsMult);
94
95 TH2F *hTRKVertexYvsMult = new TH2F("hTRKVertexYvsMult", "TRKVertex Y vs mult", 200, -1, 1, 300, 0, 3000);
96 fOutput->Add(hTRKVertexYvsMult);
bbb57828 97
98 TH2F *hTRKVertexXZ = new TH2F("hTRKVertexXZ", "TRKVertex XZ corr", 200, -1, 1, 200, -20, 20);
99 fOutput->Add(hTRKVertexXZ);
791a58dc 100
bbb57828 101 TH2F *hTRKVertexYZ = new TH2F("hTRKVertexYZ", "TRKVertex YZ corr", 200, -1, 1, 200, -20, 20);
102 fOutput->Add(hTRKVertexYZ);
103
104 TH1F* hTRKVertexXdefMult = new TH1F("hTRKVertexXdefMult","TRKVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
105 fOutput->Add(hTRKVertexXdefMult);
106 TH1F* hTRKVertexYdefMult = new TH1F("hTRKVertexYdefMult","TRKVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
107 fOutput->Add(hTRKVertexYdefMult);
a680415a 108
791a58dc 109 TH1F* hTRKVertexXHighMult = new TH1F("hTRKVertexXHighMult","TRKVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
110 fOutput->Add(hTRKVertexXHighMult);
111 TH1F* hTRKVertexYHighMult = new TH1F("hTRKVertexYHighMult","TRKVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
112 fOutput->Add(hTRKVertexYHighMult);
113
bbb57828 114 TH1F* hITSSAVertexX = new TH1F("hITSSAVertexX","ITSSAVertex x; x vertex [cm]; events",200,-1,1);
115 fOutput->Add(hITSSAVertexX);
116 TH1F* hITSSAVertexY = new TH1F("hITSSAVertexY","ITSSAVertex y; y vertex [cm]; events",200,-1,1);
117 fOutput->Add(hITSSAVertexY);
118 TH1F* hITSSAVertexZ = new TH1F("hITSSAVertexZ","ITSSAVertex z; z vertex [cm]; events",200,-20,20);
119 fOutput->Add(hITSSAVertexZ);
120
121 TH2F *hITSSAVertexXZ = new TH2F("hITSSAVertexXZ", "ITSSAVertex XZ corr", 200, -1, 1, 200, -20, 20);
122 fOutput->Add(hITSSAVertexXZ);
123
124 TH2F *hITSSAVertexYZ = new TH2F("hITSSAVertexYZ", "ITSSAVertex YZ corr", 200, -1, 1, 200, -20, 20);
125 fOutput->Add(hITSSAVertexYZ);
791a58dc 126
127 TH2F *hITSSAVertexXvsMult = new TH2F("hITSSAVertexXvsMult", "ITSSAVertex X vs mult", 200, -1, 1, 300, 0, 3000);
128 fOutput->Add(hITSSAVertexXvsMult);
129
130 TH2F *hITSSAVertexYvsMult = new TH2F("hITSSAVertexYvsMult", "ITSSAVertex Y vs mult", 200, -1, 1, 300, 0, 3000);
131 fOutput->Add(hITSSAVertexYvsMult);
bbb57828 132
133 TH1F* hITSSAVertexXdefMult = new TH1F("hITSSAVertexXdefMult","ITSSAVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
134 fOutput->Add(hITSSAVertexXdefMult);
135 TH1F* hITSSAVertexYdefMult = new TH1F("hITSSAVertexYdefMult","ITSSAVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
136 fOutput->Add(hITSSAVertexYdefMult);
137
791a58dc 138
139 TH1F* hITSSAVertexXHighMult = new TH1F("hITSSAVertexXHighMult","ITSSAVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
140 fOutput->Add(hITSSAVertexXHighMult);
141 TH1F* hITSSAVertexYHighMult = new TH1F("hITSSAVertexYHighMult","ITSSAVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
142 fOutput->Add(hITSSAVertexYHighMult);
143
bbb57828 144 PostData(1, fOutput);
145
146 return;
147}
148
149
150//________________________________________________________________________
151void AliMeanVertexCalibTask::UserExec(Option_t *)
152{
153 // Main loop
154 // Called for each event
155
156 if (!InputEvent()) {
157 Printf("ERROR: fESD not available");
158 return;
159 }
160
161 AliESDEvent* esdE = (AliESDEvent*) InputEvent();
162
bbb57828 163 const AliMultiplicity *alimult = esdE->GetMultiplicity();
164 Int_t ntrklets=0;
165 if(alimult) ntrklets = alimult->GetNumberOfTracklets();
29be8a66 166
167 const char* beamType = esdE->GetBeamType();
87b471fa 168 // Printf("beam type = %s", beamType);
29be8a66 169
170 Bool_t kLowFlux = kTRUE, kHighFlux = kFALSE;
171 // TString pp= "p-p";
172 //TString pA= "p-A";
173 TString AA= "A-A";
174
175 if (beamType == AA){
176 kHighFlux = kTRUE;
177 kLowFlux = kFALSE;
87b471fa 178 // Printf ("high flux setting");
29be8a66 179 }
180
181 AliCDBManager* man = AliCDBManager::Instance();
9fa1a9e8 182 //man->SetDefaultStorage("raw://");
29be8a66 183 Int_t runNb = esdE->GetRunNumber();
791a58dc 184 if (runNb > 0) {
185 man->SetRun(runNb);
87b471fa 186 // Printf("runNb = %d", runNb);
791a58dc 187 }
188
29be8a66 189 AliCDBEntry *entry = (AliCDBEntry*)man->Get("GRP/Calib/RecoParam/");
87b471fa 190 // Printf("entry = %p", entry);
791a58dc 191 TObjArray *arrayRecoParam=0x0;
192 if (entry) {
193 arrayRecoParam = (TObjArray*)entry->GetObject();
87b471fa 194 // Printf("arrayRecoParam = %p", arrayRecoParam);
791a58dc 195 }
196 else {
197 Printf("CDBEntry not found");
198 return;
199 }
29be8a66 200 AliGRPRecoParam *grpRecoParam=0x0;
201 if (kLowFlux) grpRecoParam= (AliGRPRecoParam*)arrayRecoParam->At(1);
202 else if (kHighFlux) grpRecoParam = (AliGRPRecoParam*)arrayRecoParam->At(2);
203
29be8a66 204 AliVertexerTracks *vertexer= new AliVertexerTracks(esdE->GetMagneticField());
205 vertexer->SetITSMode();
206 vertexer->SetConstraintOff();
207
208 if (grpRecoParam) {
209 Int_t nCutsVertexer = grpRecoParam->GetVertexerTracksNCuts();
210 Double_t *cutsVertexer = new Double_t[nCutsVertexer];
211 grpRecoParam->GetVertexerTracksCutsITS(cutsVertexer,nCutsVertexer);
212 vertexer->SetCuts(cutsVertexer,nCutsVertexer);
213 delete [] cutsVertexer; cutsVertexer = NULL;
214 }
bbb57828 215
29be8a66 216 vertexer->SetConstraintOff();
bbb57828 217
29be8a66 218 //track vertex recomputed from the vertexer
219 AliESDVertex *trkv = vertexer->FindPrimaryVertex(esdE);
87b471fa 220
29be8a66 221 //const AliESDVertex *trkv = esdE->GetPrimaryVertexTracks();
222
223 //SPD vertex taken from the ESD
bbb57828 224 const AliESDVertex *spdv = esdE->GetPrimaryVertexSPD();
225
29be8a66 226 //ITSSA vertex computed from the vertexer
227 vertexer->SetITSpureSA();
228 AliESDVertex *itsSAv = vertexer->FindPrimaryVertex(esdE);
bbb57828 229
230 if(spdv) {
231 if(spdv->GetNContributors()>0) {
232 TString title=spdv->GetTitle();
233 if(title.Contains("3D")) {
234 ((TH1F*)fOutput->FindObject("hSPDVertexX"))->Fill(spdv->GetXv());
235 ((TH1F*)fOutput->FindObject("hSPDVertexY"))->Fill(spdv->GetYv());
236 }
237 ((TH1F*)fOutput->FindObject("hSPDVertexZ"))->Fill(spdv->GetZv());
238 }
239 }
240
241
242 if(trkv) {
243 if(trkv->GetNContributors()>0) {
244 ((TH1F*)fOutput->FindObject("hTRKVertexX"))->Fill(trkv->GetXv());
245 ((TH1F*)fOutput->FindObject("hTRKVertexY"))->Fill(trkv->GetYv());
246 ((TH1F*)fOutput->FindObject("hTRKVertexZ"))->Fill(trkv->GetZv());
791a58dc 247
248 ((TH2F*)fOutput->FindObject("hTRKVertexXvsMult"))->Fill(trkv->GetXv(), ntrklets);
249 ((TH2F*)fOutput->FindObject("hTRKVertexYvsMult"))->Fill(trkv->GetYv(), ntrklets);
bbb57828 250
251 if (ntrklets>30 && ntrklets<45){
252 ((TH1F*)fOutput->FindObject("hTRKVertexXdefMult"))->Fill(trkv->GetXv());
253 ((TH1F*)fOutput->FindObject("hTRKVertexYdefMult"))->Fill(trkv->GetYv());
254 }
29be8a66 255
256 if (ntrklets>1500){
257 ((TH1F*)fOutput->FindObject("hTRKVertexXHighMult"))->Fill(trkv->GetXv());
258 ((TH1F*)fOutput->FindObject("hTRKVertexYHighMult"))->Fill(trkv->GetYv());
259 }
bbb57828 260
261 ((TH2F*)fOutput->FindObject("hTRKVertexXZ"))->Fill(trkv->GetXv(),trkv->GetZv());
262 ((TH2F*)fOutput->FindObject("hTRKVertexYZ"))->Fill(trkv->GetYv(),trkv->GetZv());
263
264 }
265 }
29be8a66 266
bbb57828 267 if (itsSAv){
268 if (itsSAv->GetNContributors()>0){
29be8a66 269
bbb57828 270 ((TH1F*)fOutput->FindObject("hITSSAVertexX"))->Fill(itsSAv->GetXv());
271 ((TH1F*)fOutput->FindObject("hITSSAVertexY"))->Fill(itsSAv->GetYv());
272 ((TH1F*)fOutput->FindObject("hITSSAVertexZ"))->Fill(itsSAv->GetZv());
791a58dc 273
274 ((TH2F*)fOutput->FindObject("hITSSAVertexXvsMult"))->Fill(itsSAv->GetXv(), ntrklets);
275 ((TH2F*)fOutput->FindObject("hITSSAVertexYvsMult"))->Fill(itsSAv->GetYv(), ntrklets);
bbb57828 276
29be8a66 277 if (ntrklets>30 && ntrklets<45){
278 ((TH1F*)fOutput->FindObject("hITSSAVertexXdefMult"))->Fill(itsSAv->GetXv());
279 ((TH1F*)fOutput->FindObject("hITSSAVertexYdefMult"))->Fill(itsSAv->GetYv());
280 }
281
282 if (ntrklets>1500){
283 ((TH1F*)fOutput->FindObject("hITSSAVertexXHighMult"))->Fill(itsSAv->GetXv());
284 ((TH1F*)fOutput->FindObject("hITSSAVertexYHighMult"))->Fill(itsSAv->GetYv());
285 }
bbb57828 286
287 ((TH2F*)fOutput->FindObject("hITSSAVertexXZ"))->Fill(itsSAv->GetXv(),itsSAv->GetZv());
288 ((TH2F*)fOutput->FindObject("hITSSAVertexYZ"))->Fill(itsSAv->GetYv(),itsSAv->GetZv());
29be8a66 289
bbb57828 290 }
291 }
292
293 delete itsSAv;
29be8a66 294 delete vertexer;
bbb57828 295
296 PostData(1, fOutput);
297
298 return;
299
300}
301
302
303//________________________________________________________________________
304void AliMeanVertexCalibTask::Terminate(Option_t *)
305{
306 // Draw result to the screen
307 // Called once at the end of the query
308 fOutput = dynamic_cast<TList*> (GetOutputData(1));
309 if (!fOutput) {
310 Printf("ERROR: fOutput not available");
311 return;
312 }
313
314
315 return;
316
317}
318
319
320//__________________________________________________________________________
29be8a66 321// AliESDVertex* AliMeanVertexCalibTask::ReconstructPrimaryVertex(Bool_t constr,Int_t mode) const {
322// // On the fly reco of ITS+TPC vertex from ESD
323// // mode=0 use all tracks
324// // mode=1 use odd-number tracks
325// // mode=2 use even-number tracks
326
327// AliESDEvent* evt = (AliESDEvent*) fInputEvent;
328// AliVertexerTracks vertexer(evt->GetMagneticField());
329// if(evt->GetNumberOfTracks()<500) {
330// vertexer.SetITSMode(); // defaults
331// vertexer.SetMinClusters(4); // default is 5
332// } else {
333// vertexer.SetITSMode(0.1,0.1,0.5,5,1,3.,100.,1000.,3.,30.,1,1);// PbPb
334// }
335// if (fOnlyITSSATracks) vertexer.SetITSpureSA();
bbb57828 336
29be8a66 337// Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
338// Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
339// Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
340// AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
341// vertexer.SetVtxStart(initVertex);
342// delete initVertex;
343// if(!constr) vertexer.SetConstraintOff();
344
345// if(fOnlyITSTPCTracks || fOnlyITSSATracks || mode>0) {
346// Int_t iskip=0;
347// Int_t ntracks = evt->GetNumberOfTracks();
348// Int_t *skip = new Int_t[ntracks];
349// for(Int_t i=0;i<ntracks;i++) skip[i]=-1;
350// for(Int_t itr=0;itr<ntracks; itr++) {
351// AliESDtrack* track = evt->GetTrack(itr);
352// if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
353// skip[iskip++]=itr;
354// continue;
355// }
356// if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
357// skip[iskip++]=itr;
358// continue;
359// }
360// if(mode==1 && itr%2==0) skip[iskip++]=itr;
361// if(mode==2 && itr%2==1) skip[iskip++]=itr;
362// }
363// vertexer.SetSkipTracks(iskip,skip);
364// delete [] skip; skip=NULL;
365// }
bbb57828 366
29be8a66 367// return vertexer.FindPrimaryVertex(evt);
368// }