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