Example jet analysis configuration.
[u/mrichter/AliRoot.git] / TPC / AliBarrelRec_TPCparam.C
CommitLineData
a654e257 1void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
056891c0 2 //
a654e257 3 // Macro to create AliESDs.root using parametrized TPC tracking
4 // and AliITStrackerV2
056891c0 5 //
a654e257 6 // Input files:
7 // - galice.root
8 // - Kinematics.root
9 // - TrackReferences.root
10 // - ITS.RecPoints.root (use AliRecontruction class)
11 // - ITS.Vertex.root (use $ALICE_ROOT/ITS/AliITSVertexerZTest.C)
e0cf7c8a 12 //
a654e257 13 // A. Dainese - LNL
e0cf7c8a 14 //
056891c0 15
4201c420 16 Int_t collcode = 1; // pp collisions
17 Bool_t useMeanVtx = kTRUE;
a654e257 18
19 if (gAlice) {
20 delete gAlice->GetRunLoader();
21 delete gAlice;
22 gAlice=0;
23 }
24
25 AliRunLoader *rl = AliRunLoader::Open("galice.root");
26 if (rl == 0x0) {
27 cerr<<"Can not open session"<<endl;
28 return;
29 }
30 Int_t retval = rl->LoadgAlice();
31 if (retval) {
32 cerr<<"LoadgAlice returned error"<<endl;
33 delete rl;
34 return;
35 }
36 retval = rl->LoadHeader();
37 if (retval) {
38 cerr<<"LoadHeader returned error"<<endl;
39 delete rl;
40 return;
41 }
42 gAlice=rl->GetAliRun();
43
44
45 TDatabasePDG *DataBase = TDatabasePDG::Instance();
46
47 // Get field from galice.root
48 AliMagF *fiel = (AliMagF*)gAlice->Field();
49 Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
50 // Set the conversion constant between curvature and Pt
51 AliTracker::SetFieldMap(fiel,kTRUE);
52
53 /**** The TPC corner ********************/
54
a654e257 55 AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
56 tpcTrackerPar.Init();
57
58 //**** Switch on the PID class (mandatory!)
59 AliPID pid;
60
61 /**** The ITS corner ********************/
62
63 AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
64 if (itsl == 0x0) {
65 cerr<<"Cannot get the ITS loader"<<endl;
66 return;
67 }
68 itsl->LoadRecPoints("read");
69
70 AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
71 if (!dITS) {
72 cerr<<"Cannot find the ITS detector !"<<endl;
73 return;
74 }
75 AliITSgeom *geom = dITS->GetITSgeom();
76
77 //An instance of the ITS tracker
78 AliITStrackerV2 itsTracker(geom);
79 Int_t ITSclusters[6] = {1,1,1,1,1,1};
80 itsTracker.SetLayersNotToSkip(ITSclusters);
81
4201c420 82 // Primary vertex reconstruction in pp
83 AliESDVertex *initVertex = 0;
84 if(collcode==1 && useMeanVtx) {
85 // open the mean vertex
86 TFile *invtx = new TFile("AliESDVertexMean.root");
87 initVertex = (AliESDVertex*)invtx->Get("vtxmean");
88 invtx->Close();
89 delete invtx;
90 } else {
91 Double_t pos[3]={0.5,0.5,0.};
92 Double_t err[3]={3.,3.,5.};
93 initVertex = new AliESDVertex(pos,err);
94 }
95 AliVertexerTracks *vertexer = new AliVertexerTracks;
96 vertexer->SetVtxStart(initVertex);
97
98
99 /***** The TREE for ESD is born *****/
a654e257 100 TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
101 AliESD *event=0;
102 esdTree->Branch("ESD","AliESD",&event);
103
104 if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
105 if(lastEvent>rl->GetNumberOfEvents()) lastEvent=rl->GetNumberOfEvents()-1;
106 cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
107
108 TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
4201c420 109 AliESDVertex *vertexSPD = new AliESDVertex();
a654e257 110 Char_t zver[100];
111 Double_t vtx[3];
112 Double_t cvtx[6];
113
114
4201c420 115 //<---------------------------------- The Loop over events begins
a654e257 116 TStopwatch timer;
117 Int_t trc;
118 for(Int_t i=firstEvent; i<=lastEvent; i++) {
119
120 cerr<<" Processing event number : "<<i<<endl;
121 AliESD *event = new AliESD();
122 event->SetRunNumber(gAlice->GetRunNumber());
123 event->SetEventNumber(i);
124 event->SetMagneticField(gAlice->Field()->SolenoidField());
125 rl->GetEvent(i);
126
4201c420 127 //***** Primary vertex from SPD from file
a654e257 128 sprintf(zver,"Event%d/Vertex",i);
4201c420 129 vertexSPD = (AliESDVertex*)ppZ->Get(zver);
130 if(!vertexSPD) {
a654e257 131 esdTree->Fill(); delete event;
132 continue;
133 }
4201c420 134 event->SetVertex(vertexSPD);
135 vertexSPD->GetXYZ(vtx);
136 vertexSPD->GetCovMatrix(cvtx);
a654e257 137
138 //***** TPC tracking
139 if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
140 printf("exiting TPC tracker with code %d in event %d\n",trc,i);
141 esdTree->Fill(); delete event;
142 continue;
143 }
144
145 //***** ITS tracking
146 itsTracker.SetVertex(vtx,cvtx);
147 TTree *itsTree=itsl->TreeR();
148 if (!itsTree) {
149 cerr<<"Can't get the ITS cluster tree !\n";
150 esdTree->Fill(); delete event;
151 return;
152 }
153 itsTracker.UnloadClusters();
154 itsTracker.LoadClusters(itsTree);
155 if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
156 printf("exiting ITS tracker with code %d in event %d\n",trc,i);
157 esdTree->Fill(); delete event;
158 continue;
159 }
160
4201c420 161 //***** Propagate kITSin tracks to local x = 0 (through beam pipe)
162 ToLocalX0(event);
163
164
165 //***** Vertex from ESD tracks
166 if(collcode==1) { // pp
167 AliESDVertex *vertexTrks =
168 (AliESDVertex*)vertexer->FindPrimaryVertex(event);
169 event->SetPrimaryVertex(vertexTrks);
170 }
a654e257 171
172 esdTree->Fill();
173 delete event;
174
175 }//<-----------------------------------The Loop over events ends here
176 timer.Stop(); timer.Print();
177
4201c420 178 // The AliESDs.root is born
a654e257 179 TFile *ef = TFile::Open("AliESDs.root","RECREATE");
180 if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
181
4201c420 182 //Write the TREE and close everything
183 esdTree->Write();
a654e257 184 delete esdTree;
185 ef->Close();
186
4201c420 187 delete vertexer;
188 delete initVertex;
a654e257 189 delete rl;
056891c0 190
e0cf7c8a 191 return;
056891c0 192}
4201c420 193//--------------------------------------------------------------------------
194Int_t ToLocalX0(AliESD *esd) {
195
196 Int_t ntracks = esd->GetNumberOfTracks();
197 AliESDtrack *esdTrack = 0;
198
199 // loop on tracks
200 for(Int_t tr = 0; tr<ntracks; tr++) {
201 esdTrack = (AliESDtrack*)esd->GetTrack(tr);
202 // ask for kITSin
203 if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
204 AliITStrackV2 *itsTrack = 0;
205 try {
206 itsTrack = new AliITStrackV2(*esdTrack);
207 esdTrack = 0;
208 }
209 catch (const Char_t *msg) {
210 Warning("FindPrimaryVertexForCurrentEvent",msg);
211 continue;
212 }
213 // propagate track to beam pipe (0.8 mm of Be)
214 if(itsTrack->GetX()>3.) itsTrack->PropagateTo(3.,0.0023,65.19);
215 // propagate track to (0,0)
216 itsTrack->PropagateTo(0.,0.,0.);
217 itsTrack->UpdateESDtrack(AliESDtrack::kITSin);
218 esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
219 delete itsTrack;
220 } // end loop on tracks
221
222 return 0;
223}