1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //*************************************************************************
17 // Class AliMeanVertexCalibTask
18 // AliAnalysisTask to extract from ESD the information on primary vertex
19 // reconstruction in order to compute the MeanVertex object
21 // Author: D.Caffarri, davide.caffarri@pd.infn.it
22 // A.Dainese, andrea.dainese@pd.infn.it
23 //*************************************************************************
30 #include "AliMultiplicity.h"
31 #include "AliESDVertex.h"
32 #include "AliESDEvent.h"
33 #include "AliVertexerTracks.h"
34 #include "AliCDBManager.h"
35 #include "AliCDBEntry.h"
36 #include "AliGRPRecoParam.h"
38 #include "AliMeanVertexCalibTask.h"
41 ClassImp(AliMeanVertexCalibTask)
43 //_____________________________________________________________________
44 AliMeanVertexCalibTask::AliMeanVertexCalibTask(const char *name):
45 AliAnalysisTaskSE(name),
48 fOnlyITSTPCTracks(kFALSE),
49 fOnlyITSSATracks(kTRUE)
54 // Define input and output slots here
55 // Output slot #0 writes into a TList container
56 DefineOutput(1, TList::Class()); //My private output
59 //________________________________________________________________________
60 AliMeanVertexCalibTask::~AliMeanVertexCalibTask()
70 //________________________________________________________________________
71 void AliMeanVertexCalibTask::UserCreateOutputObjects()
74 // Several histograms are more conveniently managed in a TList
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);
92 TH2F *hTRKVertexXvsMult = new TH2F("hTRKVertexXvsMult", "TRKVertex X vs mult", 200, -1, 1, 300, 0, 3000);
93 fOutput->Add(hTRKVertexXvsMult);
95 TH2F *hTRKVertexYvsMult = new TH2F("hTRKVertexYvsMult", "TRKVertex Y vs mult", 200, -1, 1, 300, 0, 3000);
96 fOutput->Add(hTRKVertexYvsMult);
98 TH2F *hTRKVertexXZ = new TH2F("hTRKVertexXZ", "TRKVertex XZ corr", 200, -1, 1, 200, -20, 20);
99 fOutput->Add(hTRKVertexXZ);
101 TH2F *hTRKVertexYZ = new TH2F("hTRKVertexYZ", "TRKVertex YZ corr", 200, -1, 1, 200, -20, 20);
102 fOutput->Add(hTRKVertexYZ);
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);
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);
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);
121 TH2F *hITSSAVertexXZ = new TH2F("hITSSAVertexXZ", "ITSSAVertex XZ corr", 200, -1, 1, 200, -20, 20);
122 fOutput->Add(hITSSAVertexXZ);
124 TH2F *hITSSAVertexYZ = new TH2F("hITSSAVertexYZ", "ITSSAVertex YZ corr", 200, -1, 1, 200, -20, 20);
125 fOutput->Add(hITSSAVertexYZ);
127 TH2F *hITSSAVertexXvsMult = new TH2F("hITSSAVertexXvsMult", "ITSSAVertex X vs mult", 200, -1, 1, 300, 0, 3000);
128 fOutput->Add(hITSSAVertexXvsMult);
130 TH2F *hITSSAVertexYvsMult = new TH2F("hITSSAVertexYvsMult", "ITSSAVertex Y vs mult", 200, -1, 1, 300, 0, 3000);
131 fOutput->Add(hITSSAVertexYvsMult);
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);
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);
144 PostData(1, fOutput);
150 //________________________________________________________________________
151 void AliMeanVertexCalibTask::UserExec(Option_t *)
154 // Called for each event
157 Printf("ERROR: fESD not available");
161 AliESDEvent* esdE = (AliESDEvent*) InputEvent();
163 const AliMultiplicity *alimult = esdE->GetMultiplicity();
165 if(alimult) ntrklets = alimult->GetNumberOfTracklets();
167 const char* beamType = esdE->GetBeamType();
168 // Printf("beam type = %s", beamType);
170 Bool_t kLowFlux = kTRUE, kHighFlux = kFALSE;
171 // TString pp= "p-p";
178 // Printf ("high flux setting");
181 AliCDBManager* man = AliCDBManager::Instance();
182 //man->SetDefaultStorage("raw://");
183 Int_t runNb = esdE->GetRunNumber();
186 // Printf("runNb = %d", runNb);
189 AliCDBEntry *entry = (AliCDBEntry*)man->Get("GRP/Calib/RecoParam/");
190 // Printf("entry = %p", entry);
191 TObjArray *arrayRecoParam=0x0;
193 arrayRecoParam = (TObjArray*)entry->GetObject();
194 // Printf("arrayRecoParam = %p", arrayRecoParam);
197 Printf("CDBEntry not found");
200 AliGRPRecoParam *grpRecoParam=0x0;
201 if (kLowFlux) grpRecoParam= (AliGRPRecoParam*)arrayRecoParam->At(1);
202 else if (kHighFlux) grpRecoParam = (AliGRPRecoParam*)arrayRecoParam->At(2);
204 AliVertexerTracks *vertexer= new AliVertexerTracks(esdE->GetMagneticField());
205 vertexer->SetITSMode();
206 vertexer->SetConstraintOff();
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;
216 vertexer->SetConstraintOff();
218 //track vertex recomputed from the vertexer
219 AliESDVertex *trkv = vertexer->FindPrimaryVertex(esdE);
221 //const AliESDVertex *trkv = esdE->GetPrimaryVertexTracks();
223 //SPD vertex taken from the ESD
224 const AliESDVertex *spdv = esdE->GetPrimaryVertexSPD();
226 //ITSSA vertex computed from the vertexer
227 vertexer->SetITSpureSA();
228 AliESDVertex *itsSAv = vertexer->FindPrimaryVertex(esdE);
231 if(spdv->GetNContributors()>0) {
232 TString title=spdv->GetTitle();
233 if(title.Contains("3D")) {
234 ((TH1F*)fOutput->FindObject("hSPDVertexX"))->Fill(spdv->GetX());
235 ((TH1F*)fOutput->FindObject("hSPDVertexY"))->Fill(spdv->GetY());
237 ((TH1F*)fOutput->FindObject("hSPDVertexZ"))->Fill(spdv->GetZ());
243 if(trkv->GetNContributors()>0) {
244 ((TH1F*)fOutput->FindObject("hTRKVertexX"))->Fill(trkv->GetX());
245 ((TH1F*)fOutput->FindObject("hTRKVertexY"))->Fill(trkv->GetY());
246 ((TH1F*)fOutput->FindObject("hTRKVertexZ"))->Fill(trkv->GetZ());
248 ((TH2F*)fOutput->FindObject("hTRKVertexXvsMult"))->Fill(trkv->GetX(), ntrklets);
249 ((TH2F*)fOutput->FindObject("hTRKVertexYvsMult"))->Fill(trkv->GetY(), ntrklets);
251 if (ntrklets>30 && ntrklets<45){
252 ((TH1F*)fOutput->FindObject("hTRKVertexXdefMult"))->Fill(trkv->GetX());
253 ((TH1F*)fOutput->FindObject("hTRKVertexYdefMult"))->Fill(trkv->GetY());
257 ((TH1F*)fOutput->FindObject("hTRKVertexXHighMult"))->Fill(trkv->GetX());
258 ((TH1F*)fOutput->FindObject("hTRKVertexYHighMult"))->Fill(trkv->GetY());
261 ((TH2F*)fOutput->FindObject("hTRKVertexXZ"))->Fill(trkv->GetX(),trkv->GetZ());
262 ((TH2F*)fOutput->FindObject("hTRKVertexYZ"))->Fill(trkv->GetY(),trkv->GetZ());
268 if (itsSAv->GetNContributors()>0){
270 ((TH1F*)fOutput->FindObject("hITSSAVertexX"))->Fill(itsSAv->GetX());
271 ((TH1F*)fOutput->FindObject("hITSSAVertexY"))->Fill(itsSAv->GetY());
272 ((TH1F*)fOutput->FindObject("hITSSAVertexZ"))->Fill(itsSAv->GetZ());
274 ((TH2F*)fOutput->FindObject("hITSSAVertexXvsMult"))->Fill(itsSAv->GetX(), ntrklets);
275 ((TH2F*)fOutput->FindObject("hITSSAVertexYvsMult"))->Fill(itsSAv->GetY(), ntrklets);
277 if (ntrklets>30 && ntrklets<45){
278 ((TH1F*)fOutput->FindObject("hITSSAVertexXdefMult"))->Fill(itsSAv->GetX());
279 ((TH1F*)fOutput->FindObject("hITSSAVertexYdefMult"))->Fill(itsSAv->GetY());
283 ((TH1F*)fOutput->FindObject("hITSSAVertexXHighMult"))->Fill(itsSAv->GetX());
284 ((TH1F*)fOutput->FindObject("hITSSAVertexYHighMult"))->Fill(itsSAv->GetY());
287 ((TH2F*)fOutput->FindObject("hITSSAVertexXZ"))->Fill(itsSAv->GetX(),itsSAv->GetZ());
288 ((TH2F*)fOutput->FindObject("hITSSAVertexYZ"))->Fill(itsSAv->GetY(),itsSAv->GetZ());
296 PostData(1, fOutput);
303 //________________________________________________________________________
304 void AliMeanVertexCalibTask::Terminate(Option_t *)
306 // Draw result to the screen
307 // Called once at the end of the query
308 fOutput = dynamic_cast<TList*> (GetOutputData(1));
310 Printf("ERROR: fOutput not available");
320 //__________________________________________________________________________
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
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
333 // vertexer.SetITSMode(0.1,0.1,0.5,5,1,3.,100.,1000.,3.,30.,1,1);// PbPb
335 // if (fOnlyITSSATracks) vertexer.SetITSpureSA();
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();
345 // if(fOnlyITSTPCTracks || fOnlyITSSATracks || mode>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;
356 // if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
357 // skip[iskip++]=itr;
360 // if(mode==1 && itr%2==0) skip[iskip++]=itr;
361 // if(mode==2 && itr%2==1) skip[iskip++]=itr;
363 // vertexer.SetSkipTracks(iskip,skip);
364 // delete [] skip; skip=NULL;
367 // return vertexer.FindPrimaryVertex(evt);