Merge branch 'master' into TPCdev
[u/mrichter/AliRoot.git] / ITS / DoVerticesSPD.C
CommitLineData
f5826889 1#if !defined(__CINT__) || defined(__MAKECINT__)
2#include "TFile.h"
3#include "TGeoManager.h"
4#include "TInterpreter.h"
5#include "TClassTable.h"
6#include "TNtuple.h"
7#include "TTree.h"
8#include "TArrayF.h"
9#include "TParticle.h"
10#include "AliESD.h"
11#include "AliGenPythiaEventHeader.h"
12#include "AliTracker.h"
13#include "AliHeader.h"
14#include "AliITSLoader.h"
4c00c36c 15#include "AliITSsegmentationSPD.h"
f5826889 16#include "AliVertexerTracks.h"
17#include "AliCDBManager.h"
18#include "AliGeomManager.h"
332910f8 19#include "AliGRPManager.h"
fdc48e5f 20#include "AliITSDetTypeRec.h"
f5826889 21#include "AliITSVertexer3D.h"
22#include "AliITSVertexerZ.h"
23#include "AliESDVertex.h"
7d766abb 24#include "AliESDEvent.h"
f5826889 25#include "AliRun.h"
26#include "AliRunLoader.h"
27#include "AliGenEventHeader.h"
28#include "AliStack.h"
29#endif
30
31/* $Id$ */
32
1ff24d0a 33Bool_t DoVerticesSPD(Bool_t isMC=kFALSE, Int_t pileupalgo=1, Int_t optdebug=0){
792b611a 34
f5826889 35 TFile *fint = new TFile("VertexSPD.root","recreate");
8bff6fba 36 TNtuple *nt = new TNtuple("ntvert","vertices","xtrue:ytrue:ztrue:zZ:zdiffZ:zerrZ:ntrksZ:x3D:xdiff3D:xerr3D:y3D:ydiff3D:yerr3D:z3D:zdiff3D:zerr3D:ntrks3D:nchgen:ntrklets:nrp1:ptyp:is3D:isTriggered");
f5826889 37
38 if (gClassTable->GetID("AliRun") < 0) {
39 gInterpreter->ExecuteMacro("loadlibs.C");
40 }
41 // Set OCDB if needed
42 AliCDBManager* man = AliCDBManager::Instance();
43 if (!man->IsDefaultStorageSet()) {
44 printf("Setting a local default storage and run number 0\n");
162637e4 45 man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
4c00c36c 46 // man->SetSpecificStorage("GRP/GRP/Data",Form("local://%s",gSystem->pwd()));
f5826889 47 man->SetRun(0);
48 }
49 else {
50 printf("Using deafult storage \n");
51 }
52
53 // retrives geometry
54 if(!gGeoManager){
55 AliGeomManager::LoadGeometry("geometry.root");
56 }
3714ff25 57 AliGeomManager::ApplyAlignObjsFromCDB("ITS");
f5826889 58
59 AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
60 if (!runLoader) {
61 Error("DoVertices", "getting run loader from file %s failed",
62 "galice.root");
63 return kFALSE;
64 }
1ff24d0a 65
66 if(isMC){
67 runLoader->LoadgAlice();
68 gAlice = runLoader->GetAliRun();
69 if (!gAlice) {
70 Error("DoVertices", "no galice object found");
71 return kFALSE;
72 }
73 runLoader->LoadKinematics();
74 runLoader->LoadHeader();
f5826889 75 }
f5826889 76 AliITSLoader* ITSloader = (AliITSLoader*) runLoader->GetLoader("ITSLoader");
77 ITSloader->LoadRecPoints("read");
78
3714ff25 79 Int_t totev=runLoader->GetNumberOfEvents();
f5826889 80 if(optdebug) printf("Number of events= %d\n",totev);
81
5951029f 82 TFile* esdFile = TFile::Open("AliESDs.root");
83 if (!esdFile || !esdFile->IsOpen()) {
84 printf("Error in opening ESD file");
85 return kFALSE;
86 }
87
88 AliESDEvent * esd = new AliESDEvent;
89 TTree* tree = (TTree*) esdFile->Get("esdTree");
90 if (!tree) {
91 printf("Error: no ESD tree found");
92 return kFALSE;
93 }
94 esd->ReadFromTree(tree);
332910f8 95 AliGRPManager * grpMan=new AliGRPManager();
96 grpMan->ReadGRPEntry();
97 grpMan->SetMagField();
98 printf("Magnetic field set to %f T\n",AliTracker::GetBz()/10.);
5951029f 99
100 AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
4c00c36c 101 AliITSsegmentation* seg = new AliITSsegmentationSPD();
102 detTypeRec->SetSegmentationModel(0,seg);
103
f5826889 104 Double_t xnom=0.,ynom=0.;
308c2f7c 105 AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
106 vertz->Init("default");
1ff24d0a 107 AliITSVertexer3D *vert3d = new AliITSVertexer3D();
308c2f7c 108 vert3d->Init("default");
792b611a 109 vert3d->SetWideFiducialRegion(40.,1.);
7d766abb 110 vert3d->SetPileupAlgo(pileupalgo);
792b611a 111 vert3d->PrintStatus();
fdc48e5f 112 vertz->SetDetTypeRec(detTypeRec);
113 vert3d->SetDetTypeRec(detTypeRec);
1ff24d0a 114 vert3d->SetComputeMultiplicity(kTRUE);
ecce5370 115 /* uncomment these lines to use diamond constrain */
116// Double_t posdiam[3]={0.03,0.1,0.};
117// Double_t sigdiam[3]={0.01,0.01,10.0};
118// AliESDVertex* diam=new AliESDVertex(posdiam,sigdiam);
119// vertz->SetVtxStart(diam);
120// vert3d->SetVtxStart(diam);
121 /* end lines to be uncommented to use diamond constrain */
f5826889 122
123 Int_t goodz=0,good3d=0;
7d766abb 124 // Trigger mask
125 ULong64_t spdFO = (1 << 14);
126 ULong64_t v0left = (1 << 11);
127 ULong64_t v0right = (1 << 12);
f5826889 128
129 for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
130 TArrayF mcVertex(3);
131 runLoader->GetEvent(iEvent);
8bff6fba 132 Double_t nChGen = 0.;
1ff24d0a 133 Int_t ptype = 0;
f5826889 134 if(optdebug){
135 printf("==============================================================\n");
1ff24d0a 136 printf("Event: %d \n",iEvent);
f5826889 137 }
1ff24d0a 138 if(isMC){
139 runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
140 AliGenPythiaEventHeader *evh=(AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader();
141 ptype = evh->ProcessType();
142
143 AliStack* stack = runLoader->Stack();
144 TTree *treek=(TTree*)runLoader->TreeK();
145 Int_t npart = (Int_t)treek->GetEntries();
146 if(optdebug) printf("Process Type = %d --- Particles %d\n",ptype,npart);
147
148 // loop on particles to get generated dN/dy
149 for(Int_t iPart=0; iPart<npart; iPart++) {
150 if(!stack->IsPhysicalPrimary(iPart)) continue;
151 TParticle* part = (TParticle*)stack->Particle(iPart);
152 if(part->GetPDG()->Charge() == 0) continue;
153 Double_t eta=part->Eta();
154
8bff6fba 155 if(TMath::Abs(eta)<1.5) nChGen+=1.;
1ff24d0a 156 }
8bff6fba 157 if(optdebug) printf(" N. of generated charged primaries in |eta|<1.5 = %f\n",nChGen);
f5826889 158 }
7d766abb 159 tree->GetEvent(iEvent);
308c2f7c 160
7d766abb 161 TTree* cltree = ITSloader->TreeR();
162 ULong64_t triggerMask=esd->GetTriggerMask();
163 // MB1: SPDFO || V0L || V0R
164 Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)));
308c2f7c 165 AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
ecce5370 166 AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);
ecce5370 167 AliMultiplicity *alimult = vert3d->GetMultiplicity();
792b611a 168 Int_t ntrklets=0;
169 Int_t nrecp1=0;
170 if(alimult){
171 ntrklets=alimult->GetNumberOfTracklets() ;
172 nrecp1=ntrklets+alimult->GetNumberOfSingleClusters();
f5826889 173 }
ecce5370 174
f5826889 175
176 TDirectory *current = gDirectory;
177 fint->cd();
7d766abb 178 Float_t xnt[23];
f5826889 179
180 Double_t zz = 0.;
181 Double_t zresz = 0.;
182 Double_t zdiffz = 0.;
183
184 Int_t ntrkz = 0;
185 if(vtxz){
186
187 ntrkz = vtxz->GetNContributors();
188 if(ntrkz>0)goodz++;
189 zz=vtxz->GetZv();
190 zresz =vtxz->GetZRes(); // microns
191 zdiffz = 10000.*(zz-mcVertex[2]); // microns
192 }
193
194 Double_t x3d = 0.;
195 Double_t xerr3d = 0.;
196 Double_t xdiff3d = 0.;
197
198 Double_t y3d = 0.;
199 Double_t yerr3d = 0.;
200 Double_t ydiff3d = 0.;
201
202 Double_t z3d = 0.;
203 Double_t zerr3d = 0.;
204 Double_t zdiff3d = 0.;
205
206 Int_t ntrk3d = 0;
6b4d9537 207 Bool_t is3d=kFALSE;
f5826889 208 if(vtx3d){
7d766abb 209
6b4d9537 210 if(vtx3d->IsFromVertexer3D()) is3d=kTRUE;
f5826889 211 ntrk3d = vtx3d->GetNContributors();
792b611a 212 if(is3d && ntrk3d>0)good3d++;
f5826889 213 x3d = vtx3d->GetXv();
214 xerr3d=vtx3d->GetXRes();
215 xdiff3d = 10000.*(x3d-mcVertex[0]); // microns
216
217 y3d = vtx3d->GetYv();
218 yerr3d=vtx3d->GetYRes();
219 ydiff3d = 10000.*(y3d-mcVertex[1]); // microns
220
221 z3d = vtx3d->GetZv();
222 zerr3d=vtx3d->GetZRes();
223 zdiff3d = 10000.*(z3d-mcVertex[2]); // microns
224
225 }
226
227 xnt[0]=mcVertex[0];//x
228 xnt[1]=mcVertex[1];//y
229 xnt[2]=mcVertex[2];//z
230
231 xnt[3]=zz;
232 xnt[4]=zdiffz;
233 xnt[5]=zresz;
234 xnt[6]=ntrkz;
235
236 xnt[7]=x3d;
237 xnt[8]=xdiff3d;
238 xnt[9]=xerr3d;
8686a665 239 xnt[10]=y3d;
f5826889 240 xnt[11]=ydiff3d;
241 xnt[12]=yerr3d;
242 xnt[13]=z3d;
243 xnt[14]=zdiff3d;
244 xnt[15]=zerr3d;
245 xnt[16]=ntrk3d;
246
8bff6fba 247 xnt[17]=nChGen;
f5826889 248 xnt[18]=float(ntrklets);
249 xnt[19]=float(nrecp1);
250 xnt[20]=float(ptype);
792b611a 251 xnt[21]=float(is3d);
7d766abb 252 xnt[22]=float(eventTriggered);
f5826889 253 nt->Fill(xnt);
254 current->cd();
255
256 if(optdebug){
1ff24d0a 257 printf("Vertexer3D: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tContributors=%d \n",x3d,y3d,z3d,ntrk3d);
258 printf("VertexerZ: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tContributors=%d \n",0.,0.,zz,ntrkz);
8bff6fba 259 if(isMC) printf("True Pos.: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tdN/dy=%.1f\n",mcVertex[0],mcVertex[1],mcVertex[2],nChGen);
1ff24d0a 260 printf("Multiplicity: Tracklets=%d ClustersLay1=%d\n",ntrklets,nrecp1);
f5826889 261 }
262
263 }
264 fint->cd();
265 nt->Write();
266 fint->Close();
267 delete fint;
268 if(optdebug){
269 printf("***********************************************\n");
270 printf("Number of Z vertices found= %d\n",goodz);
271 printf("efficiency=%f\n",float(goodz)/float(totev));
272 printf("Number of 3D vertices found= %d\n",good3d);
273 printf("efficiency=%f\n",float(good3d)/float(totev));
274 }
275 return kTRUE;
276}