aee815b5e829597011d9f4fa0400f49206e57c9e
[u/mrichter/AliRoot.git] / ITS / DoVerticesSPD.C
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"
18 #include "AliITSDetTypeRec.h"
19 #include "AliITSVertexer3D.h"
20 #include "AliITSVertexerZ.h"
21 #include "AliESDVertex.h"
22 #include "AliRun.h"
23 #include "AliRunLoader.h"
24 #include "AliGenEventHeader.h"
25 #include "AliStack.h"
26 #endif
27
28 /*  $Id$    */
29
30 Bool_t DoVerticesSPD(Int_t optdebug=1){
31
32   TFile *fint = new TFile("VertexSPD.root","recreate");
33   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");
34
35   if (gClassTable->GetID("AliRun") < 0) {
36     gInterpreter->ExecuteMacro("loadlibs.C");
37   }
38   // Set OCDB if needed
39   AliCDBManager* man = AliCDBManager::Instance();
40   if (!man->IsDefaultStorageSet()) {
41     printf("Setting a local default storage and run number 0\n");
42     man->SetDefaultStorage("local://$ALICE_ROOT");
43     man->SetRun(0);
44   }
45   else {
46     printf("Using deafult storage \n");
47   }
48  
49   // retrives geometry 
50   if(!gGeoManager){
51     AliGeomManager::LoadGeometry("geometry.root");
52   }
53
54   AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
55   if (!runLoader) {
56     Error("DoVertices", "getting run loader from file %s failed", 
57           "galice.root");
58     return kFALSE;
59   }
60   runLoader->LoadgAlice();
61   gAlice = runLoader->GetAliRun();
62   if (!gAlice) {
63     Error("DoVertices", "no galice object found");
64     return kFALSE;
65   }
66   runLoader->LoadKinematics();
67   runLoader->LoadHeader();
68   AliITSLoader* ITSloader =  (AliITSLoader*) runLoader->GetLoader("ITSLoader");
69   ITSloader->LoadRecPoints("read");
70
71   AliMagF * magf = gAlice->Field();
72   AliTracker::SetFieldMap(magf,kTRUE);
73   if(optdebug) printf("MagneticField=%f\n",AliTracker::GetBz());
74   
75   Int_t totev=runLoader->GetNumberOfEvents();
76   if(optdebug)  printf("Number of events= %d\n",totev);
77
78   AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
79  
80   Double_t xnom=0.,ynom=0.;
81   AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
82   vertz->Init("default");
83   AliITSVertexer3D *vert3d = new AliITSVertexer3D();
84   vert3d->Init("default");
85   vert3d->SetWideFiducialRegion(40.,1.);
86   vert3d->PrintStatus();
87   vertz->SetDetTypeRec(detTypeRec);
88   vert3d->SetDetTypeRec(detTypeRec);
89
90   /* uncomment these lines to use diamond constrain */
91 //   Double_t posdiam[3]={0.03,0.1,0.};
92 //   Double_t sigdiam[3]={0.01,0.01,10.0};
93 //   AliESDVertex* diam=new AliESDVertex(posdiam,sigdiam);
94 //   vertz->SetVtxStart(diam);
95 //   vert3d->SetVtxStart(diam);
96   /* end lines to be uncommented to use diamond constrain */
97
98   Int_t goodz=0,good3d=0;
99
100   for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
101     TArrayF mcVertex(3); 
102     runLoader->GetEvent(iEvent);
103     runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
104     AliGenPythiaEventHeader *evh=(AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader();
105     Int_t ptype = evh->ProcessType();
106     if(optdebug){
107       printf("==============================================================\n");
108       printf("\nEvent: %d ---- Process Type = %d \n",iEvent,ptype);
109     }
110
111     AliStack* stack = runLoader->Stack();
112     TTree *treek=(TTree*)runLoader->TreeK();
113     Int_t npart = (Int_t)treek->GetEntries();
114     if(optdebug) printf("particles  %d\n",npart);
115
116     Double_t dNchdy = 0.;
117
118    // loop on particles to get generated dN/dy
119     for(Int_t iPart=0; iPart<npart; iPart++) {
120       if(!stack->IsPhysicalPrimary(iPart)) continue;
121       TParticle* part = (TParticle*)stack->Particle(iPart);
122       if(part->GetPDG()->Charge() == 0) continue;
123       Double_t eta=part->Eta();
124
125       if(TMath::Abs(eta)<1.5) dNchdy+=1.; 
126     }
127     if(optdebug) printf(" dNch/dy = %f\n",dNchdy);
128  
129     TTree* cltree = ITSloader->TreeR();
130
131     AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
132     AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);
133     AliMultiplicity *alimult = vert3d->GetMultiplicity();
134     Int_t  ntrklets=0;
135     Int_t nrecp1=0;
136     if(alimult){
137       ntrklets=alimult->GetNumberOfTracklets() ;
138       nrecp1=ntrklets+alimult->GetNumberOfSingleClusters();
139     }
140     TString tit=vtx3d->GetTitle();
141     Bool_t is3d=kFALSE;
142     if(tit.Contains("3D")) is3d=kTRUE;
143        
144
145
146     TDirectory *current = gDirectory;
147     fint->cd();
148     Float_t xnt[22];
149
150     Double_t zz = 0.;
151     Double_t zresz = 0.;
152     Double_t zdiffz = 0.;
153
154     Int_t ntrkz = 0;
155     if(vtxz){
156  
157       ntrkz = vtxz->GetNContributors();
158       if(ntrkz>0)goodz++;
159       zz=vtxz->GetZv();
160       zresz =vtxz->GetZRes(); // microns
161       zdiffz = 10000.*(zz-mcVertex[2]); // microns
162     }
163
164     Double_t x3d = 0.;
165     Double_t xerr3d = 0.;
166     Double_t xdiff3d = 0.;
167
168     Double_t y3d = 0.;
169     Double_t yerr3d = 0.;
170     Double_t ydiff3d = 0.;
171
172     Double_t z3d = 0.;
173     Double_t zerr3d = 0.;
174     Double_t zdiff3d = 0.;
175
176     Int_t ntrk3d = 0;
177     if(vtx3d){
178
179       ntrk3d = vtx3d->GetNContributors();
180       if(is3d && ntrk3d>0)good3d++;
181       x3d = vtx3d->GetXv();
182       xerr3d=vtx3d->GetXRes();
183       xdiff3d = 10000.*(x3d-mcVertex[0]);  // microns
184
185       y3d = vtx3d->GetYv();
186       yerr3d=vtx3d->GetYRes();
187       ydiff3d = 10000.*(y3d-mcVertex[1]);  // microns
188
189       z3d = vtx3d->GetZv();
190       zerr3d=vtx3d->GetZRes();
191       zdiff3d = 10000.*(z3d-mcVertex[2]);  // microns
192
193     }
194
195     xnt[0]=mcVertex[0];//x
196     xnt[1]=mcVertex[1];//y
197     xnt[2]=mcVertex[2];//z
198
199     xnt[3]=zz;
200     xnt[4]=zdiffz;
201     xnt[5]=zresz;
202     xnt[6]=ntrkz;
203
204     xnt[7]=x3d;
205     xnt[8]=xdiff3d;
206     xnt[9]=xerr3d;
207     xnt[10]=y3d;
208     xnt[11]=ydiff3d;
209     xnt[12]=yerr3d;
210     xnt[13]=z3d;
211     xnt[14]=zdiff3d;
212     xnt[15]=zerr3d;
213     xnt[16]=ntrk3d;
214
215     xnt[17]=dNchdy;
216     xnt[18]=float(ntrklets);
217     xnt[19]=float(nrecp1);
218     xnt[20]=float(ptype);
219     xnt[21]=float(is3d);
220     nt->Fill(xnt);
221     current->cd();
222     
223     if(optdebug){
224       printf("\nVertexerZ:  \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",zz,ntrkz,mcVertex[2],zdiffz);
225       printf("Vertexer3D: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",z3d,ntrk3d,mcVertex[2],zdiff3d);
226     }
227     
228   }
229   fint->cd();
230   nt->Write();
231   fint->Close();
232   delete fint;
233   if(optdebug){
234    printf("***********************************************\n");
235    printf("Number of Z vertices found= %d\n",goodz);
236    printf("efficiency=%f\n",float(goodz)/float(totev));
237    printf("Number of 3D vertices found= %d\n",good3d);
238    printf("efficiency=%f\n",float(good3d)/float(totev));
239   }
240   return kTRUE;
241 }