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