Fix for FMD DA
[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 "AliITSsegmentationSPD.h"
16 #include "AliVertexerTracks.h"
17 #include "AliCDBManager.h"
18 #include "AliGeomManager.h"
19 #include "AliGRPManager.h"
20 #include "AliITSDetTypeRec.h"
21 #include "AliITSVertexer3D.h"
22 #include "AliITSVertexerZ.h"
23 #include "AliESDVertex.h"
24 #include "AliESDEvent.h"
25 #include "AliRun.h"
26 #include "AliRunLoader.h"
27 #include "AliGenEventHeader.h"
28 #include "AliStack.h"
29 #endif
30
31 /*  $Id$    */
32
33 Bool_t DoVerticesSPD(Bool_t isMC=kFALSE, Int_t pileupalgo=1, Int_t optdebug=0){
34
35   TFile *fint = new TFile("VertexSPD.root","recreate");
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");
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");
45     man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
46     //    man->SetSpecificStorage("GRP/GRP/Data",Form("local://%s",gSystem->pwd()));
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   }
57   AliGeomManager::ApplyAlignObjsFromCDB("ITS");
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   }
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();
75   }
76   AliITSLoader* ITSloader =  (AliITSLoader*) runLoader->GetLoader("ITSLoader");
77   ITSloader->LoadRecPoints("read");
78
79   Int_t totev=runLoader->GetNumberOfEvents();
80   if(optdebug)  printf("Number of events= %d\n",totev);
81
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);
95   AliGRPManager * grpMan=new AliGRPManager();
96   grpMan->ReadGRPEntry();
97   grpMan->SetMagField();
98   printf("Magnetic field set to %f T\n",AliTracker::GetBz()/10.);
99
100   AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
101   AliITSsegmentation* seg = new AliITSsegmentationSPD();  
102   detTypeRec->SetSegmentationModel(0,seg);
103
104   Double_t xnom=0.,ynom=0.;
105   AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
106   vertz->Init("default");
107   AliITSVertexer3D *vert3d = new AliITSVertexer3D();  
108   vert3d->Init("default");
109   vert3d->SetWideFiducialRegion(40.,1.);
110   vert3d->SetPileupAlgo(pileupalgo);
111   vert3d->PrintStatus();
112   vertz->SetDetTypeRec(detTypeRec);
113   vert3d->SetDetTypeRec(detTypeRec);
114   vert3d->SetComputeMultiplicity(kTRUE);
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 */
122
123   Int_t goodz=0,good3d=0;
124   // Trigger mask
125   ULong64_t spdFO = (1 << 14);
126   ULong64_t v0left = (1 << 11);
127   ULong64_t v0right = (1 << 12);
128
129   for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
130     TArrayF mcVertex(3); 
131     runLoader->GetEvent(iEvent);
132     Double_t nChGen = 0.;
133     Int_t ptype = 0;
134     if(optdebug){
135       printf("==============================================================\n");
136       printf("Event: %d \n",iEvent);
137     }
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         
155         if(TMath::Abs(eta)<1.5) nChGen+=1.; 
156       }
157       if(optdebug) printf(" N. of generated charged primaries in |eta|<1.5 = %f\n",nChGen);
158     }
159     tree->GetEvent(iEvent);
160
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)));
165     AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
166     AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);
167     AliMultiplicity *alimult = vert3d->GetMultiplicity();
168     Int_t  ntrklets=0;
169     Int_t nrecp1=0;
170     if(alimult){
171       ntrklets=alimult->GetNumberOfTracklets() ;
172       nrecp1=ntrklets+alimult->GetNumberOfSingleClusters();
173     }
174
175
176     TDirectory *current = gDirectory;
177     fint->cd();
178     Float_t xnt[23];
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->GetZ();
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;
207     Bool_t is3d=kFALSE;
208     if(vtx3d){
209
210       if(vtx3d->IsFromVertexer3D()) is3d=kTRUE;
211       ntrk3d = vtx3d->GetNContributors();
212       if(is3d && ntrk3d>0)good3d++;
213       x3d = vtx3d->GetX();
214       xerr3d=vtx3d->GetXRes();
215       xdiff3d = 10000.*(x3d-mcVertex[0]);  // microns
216
217       y3d = vtx3d->GetY();
218       yerr3d=vtx3d->GetYRes();
219       ydiff3d = 10000.*(y3d-mcVertex[1]);  // microns
220
221       z3d = vtx3d->GetZ();
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;
239     xnt[10]=y3d;
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
247     xnt[17]=nChGen;
248     xnt[18]=float(ntrklets);
249     xnt[19]=float(nrecp1);
250     xnt[20]=float(ptype);
251     xnt[21]=float(is3d);
252     xnt[22]=float(eventTriggered);
253     nt->Fill(xnt);
254     current->cd();
255     
256     if(optdebug){
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);
259       if(isMC) printf("True Pos.: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tdN/dy=%.1f\n",mcVertex[0],mcVertex[1],mcVertex[2],nChGen);
260       printf("Multiplicity: Tracklets=%d  ClustersLay1=%d\n",ntrklets,nrecp1);
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 }