Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 : /* $Id$ */
16 :
17 : /* History of cvs commits:
18 : *
19 : * $Log$
20 : * Revision 1.6 2007/08/22 09:20:50 hristov
21 : * Updated QA classes (Yves)
22 : *
23 : * Revision 1.5 2007/07/11 13:43:30 hristov
24 : * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
25 : *
26 : * Revision 1.4 2007/05/04 14:49:29 policheh
27 : * AliPHOSRecPoint inheritance from AliCluster
28 : *
29 : * Revision 1.3 2007/04/25 19:39:42 kharlov
30 : * Track extracpolation improved
31 : *
32 : * Revision 1.2 2007/04/01 19:16:52 kharlov
33 : * D.P.: Produce EMCTrackSegments using TPC/ITS tracks (no CPV)
34 : *
35 : *
36 : */
37 :
38 : //_________________________________________________________________________
39 : // Implementation version 2 of algorithm class to construct PHOS track segments
40 : // Track segment for PHOS is list of
41 : // EMC RecPoint + (possibly) projection of TPC track
42 : // To find TrackSegments we do the following:
43 : // for each EMC RecPoints we look at
44 : // TPC projections radius fRtpc.
45 : // If there is such a track
46 : // we make "Link" it is just indexes of EMC and TPC track and distance
47 : // between them in the PHOS plane.
48 : // Then we sort "Links" and starting from the
49 : // least "Link" pointing to the unassined EMC and TPC assing them to
50 : // new TrackSegment.
51 : // If there is no TPC track we make TrackSegment
52 : // consisting from EMC alone. There is no TrackSegments without EMC RecPoint.
53 : //// In principle this class should be called from AliPHOSReconstructor, but
54 : // one can use it as well in standalone mode.
55 : // Use case:
56 : // root [0] AliPHOSTrackSegmentMakerv2 * t = new AliPHOSTrackSegmentMaker("galice.root", "tracksegmentsname", "recpointsname")
57 : // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
58 : // // reads gAlice from header file "galice.root", uses recpoints stored in the branch names "recpointsname" (default = "Default")
59 : // // and saves recpoints in branch named "tracksegmentsname" (default = "recpointsname")
60 : // root [1] t->ExecuteTask()
61 : // root [3] t->SetTrackSegmentsBranch("max distance 5 cm")
62 : // root [4] t->ExecuteTask("deb all time")
63 : //
64 : //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH)
65 : //
66 :
67 : // --- ROOT system ---
68 : #include "TFile.h"
69 : #include "TTree.h"
70 : #include "TBenchmark.h"
71 :
72 : // --- Standard library ---
73 : #include "Riostream.h"
74 : // --- AliRoot header files ---
75 : #include "AliPHOSGeometry.h"
76 : #include "AliPHOSTrackSegmentMakerv2.h"
77 : #include "AliPHOSTrackSegment.h"
78 : #include "AliPHOSLink.h"
79 : #include "AliPHOSEmcRecPoint.h"
80 : #include "AliPHOSCpvRecPoint.h"
81 :
82 : #include "AliESDEvent.h"
83 : #include "AliESDtrack.h"
84 :
85 20 : ClassImp( AliPHOSTrackSegmentMakerv2)
86 :
87 :
88 : //____________________________________________________________________________
89 : AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2() :
90 0 : AliPHOSTrackSegmentMaker(),
91 0 : fDefaultInit(kTRUE),
92 0 : fWrite(kFALSE),
93 0 : fNTrackSegments(0),
94 0 : fRtpc(0.f),
95 0 : fVtx(0.,0.,0.),
96 0 : fLinkUpArray(0),
97 0 : fTPCtracks(),
98 0 : fEmcFirst(0),
99 0 : fEmcLast(0),
100 0 : fModule(0),
101 0 : fTrackSegments(NULL)
102 0 : {
103 : // default ctor (to be used mainly by Streamer)
104 :
105 0 : for(Int_t i=0; i<5; i++)fNtpcTracks[i]=0 ;
106 0 : InitParameters() ;
107 0 : }
108 :
109 : //____________________________________________________________________________
110 : AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(AliPHOSGeometry *geom) :
111 0 : AliPHOSTrackSegmentMaker(geom),
112 0 : fDefaultInit(kFALSE),
113 0 : fWrite(kFALSE),
114 0 : fNTrackSegments(0),
115 0 : fRtpc(0.f),
116 0 : fVtx(0.,0.,0.),
117 0 : fLinkUpArray(0),
118 0 : fTPCtracks(),
119 0 : fEmcFirst(0),
120 0 : fEmcLast(0),
121 0 : fModule(0),
122 0 : fTrackSegments(NULL)
123 0 : {
124 : // ctor
125 0 : for(Int_t i=0; i<5; i++)fNtpcTracks[i]=0 ;
126 :
127 0 : InitParameters() ;
128 0 : Init() ;
129 0 : fESD = 0;
130 0 : }
131 :
132 : //____________________________________________________________________________
133 : AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const AliPHOSTrackSegmentMakerv2 & tsm) :
134 0 : AliPHOSTrackSegmentMaker(tsm),
135 0 : fDefaultInit(kFALSE),
136 0 : fWrite(kFALSE),
137 0 : fNTrackSegments(0),
138 0 : fRtpc(0.f),
139 0 : fVtx(0.,0.,0.),
140 0 : fLinkUpArray(0),
141 0 : fTPCtracks(),
142 0 : fEmcFirst(0),
143 0 : fEmcLast(0),
144 0 : fModule(0),
145 0 : fTrackSegments(NULL)
146 0 : {
147 : // cpy ctor: no implementation yet
148 : // requested by the Coding Convention
149 0 : Fatal("cpy ctor", "not implemented") ;
150 0 : }
151 :
152 :
153 : //____________________________________________________________________________
154 : AliPHOSTrackSegmentMakerv2::~AliPHOSTrackSegmentMakerv2()
155 0 : {
156 : // dtor
157 : // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
158 0 : if (!fDefaultInit)
159 0 : delete fLinkUpArray ;
160 0 : if (fTrackSegments) {
161 0 : fTrackSegments->Delete();
162 0 : delete fTrackSegments;
163 : }
164 0 : }
165 :
166 : //____________________________________________________________________________
167 : void AliPHOSTrackSegmentMakerv2::FillOneModule()
168 : {
169 : // Finds first and last indexes between which
170 : // clusters from one PHOS module are
171 :
172 : //First EMC clusters
173 0 : Int_t totalEmc = fEMCRecPoints->GetEntriesFast() ;
174 0 : for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
175 0 : ((static_cast<AliPHOSRecPoint *>(fEMCRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule );
176 0 : fEmcLast ++) ;
177 :
178 : //Now TPC tracks
179 0 : if(fESD){
180 : //Do it ones, only first time
181 0 : if(fModule==1){
182 0 : Int_t nTracks = fESD->GetNumberOfTracks();
183 :
184 0 : Int_t nPHOSmod = fGeom->GetNModules() ;
185 0 : if(fTPCtracks[0].size()<(UInt_t)nTracks){
186 0 : for(Int_t i=0; i<nPHOSmod; i++)
187 0 : fTPCtracks[i].resize(nTracks) ;
188 0 : }
189 0 : for(Int_t i=0; i<5; i++)fNtpcTracks[i]=0 ;
190 0 : TVector3 inPHOS ;
191 :
192 : //In this particular case we use fixed vertex position at zero
193 0 : Double_t vtx[3]={0.,0.,0.} ;
194 : AliESDtrack *track;
195 0 : Double_t xyz[3] ;
196 0 : Double_t rEMC = fGeom->GetIPtoCrystalSurface() ; //Use here ideal geometry
197 0 : for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
198 0 : track = fESD->GetTrack(iTrack);
199 0 : for(Int_t iTestMod=1; iTestMod<=nPHOSmod; iTestMod++){
200 0 : Double_t modAngle=270.+fGeom->GetPHOSAngle(iTestMod) ;
201 0 : modAngle*=TMath::Pi()/180. ;
202 0 : track->Rotate(modAngle);
203 0 : if (!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz))
204 0 : continue; //track coord on the cylinder of PHOS radius
205 0 : if ((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0)
206 0 : continue;
207 : //Check if this track hits PHOS
208 0 : inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
209 0 : Int_t modNum ;
210 0 : Double_t x,z ;
211 0 : fGeom->ImpactOnEmc(vtx, inPHOS.Theta(), inPHOS.Phi(), modNum, z, x) ;
212 0 : if(modNum==iTestMod){
213 : //Mark this track as one belonging to module
214 0 : TrackInPHOS_t &t = fTPCtracks[modNum-1][fNtpcTracks[modNum-1]++] ;
215 0 : t.track = track ;
216 0 : t.x = x ;
217 0 : t.z = z ;
218 : break ;
219 : }
220 0 : }
221 : }
222 0 : }
223 : }
224 :
225 0 : }
226 :
227 : //____________________________________________________________________________
228 : void AliPHOSTrackSegmentMakerv2::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,
229 : AliESDtrack *track,
230 : Float_t &dx, Float_t &dz) const
231 : {
232 : // Calculates the distance between the EMC RecPoint and the CPV RecPoint
233 : // Clusters are sorted in "rows" and "columns" of width 1 cm
234 :
235 0 : TVector3 emcGlobal; // Global position of the EMC recpoint
236 0 : fGeom->GetGlobalPHOS((AliPHOSRecPoint*)emcClu,emcGlobal);
237 :
238 : //Calculate actual distance to the PHOS surface
239 0 : Double_t modAngle=270.+fGeom->GetPHOSAngle(emcClu->GetPHOSMod()) ;
240 0 : modAngle*=TMath::Pi()/180. ;
241 0 : Double_t rEMC = emcGlobal.X()*TMath::Cos(modAngle)+emcGlobal.Y()*TMath::Sin(modAngle) ;
242 0 : track->Rotate(modAngle);
243 0 : Double_t xyz[3] ;
244 0 : if(!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz)){
245 0 : dx=999. ;
246 0 : dz=999. ;
247 0 : return ; //track coord on the cylinder of PHOS radius
248 : }
249 0 : if((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0){
250 0 : dx=999. ;
251 0 : dz=999. ;
252 0 : return ;
253 : }
254 :
255 0 : dx=(emcGlobal.X()-xyz[0])*TMath::Sin(modAngle)-(emcGlobal.Y()-xyz[1])*TMath::Cos(modAngle) ;
256 0 : dx=TMath::Sign(dx,(Float_t)(emcGlobal.X()-xyz[0])) ; //set direction
257 0 : dz=emcGlobal.Z()-xyz[2] ;
258 :
259 0 : return ;
260 0 : }
261 : //____________________________________________________________________________
262 : void AliPHOSTrackSegmentMakerv2::Init()
263 : {
264 : // Make all memory allocations that are not possible in default constructor
265 :
266 0 : fLinkUpArray = new TClonesArray("AliPHOSLink", 1000);
267 0 : fTrackSegments = new TClonesArray("AliPHOSTrackSegment",100);
268 0 : }
269 :
270 : //____________________________________________________________________________
271 : void AliPHOSTrackSegmentMakerv2::InitParameters()
272 : {
273 : //Initializes parameters
274 0 : fRtpc = 4. ;
275 0 : fEmcFirst = 0 ;
276 0 : fEmcLast = 0 ;
277 0 : fLinkUpArray = 0 ;
278 0 : fWrite = kTRUE ;
279 0 : }
280 :
281 :
282 : //____________________________________________________________________________
283 : void AliPHOSTrackSegmentMakerv2::MakeLinks()
284 : {
285 : // Finds distances (links) between all EMC and CPV clusters,
286 : // which are not further apart from each other than fRcpv
287 : // and sort them in accordance with this distance
288 :
289 0 : fLinkUpArray->Clear() ;
290 :
291 : AliPHOSEmcRecPoint * emcclu ;
292 :
293 : Int_t iLinkUp = 0 ;
294 :
295 : Int_t iEmcRP;
296 0 : for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
297 0 : emcclu = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP)) ;
298 0 : TVector3 vecEmc ;
299 0 : emcclu->GetLocalPosition(vecEmc) ;
300 0 : Int_t mod=emcclu->GetPHOSMod() ;
301 0 : for(Int_t itr=0; itr<fNtpcTracks[mod-1] ; itr++){
302 0 : TrackInPHOS_t &t = fTPCtracks[mod-1][itr] ;
303 : //calculate raw distance
304 0 : Double_t rawdx = t.x - vecEmc.X() ;
305 0 : Double_t rawdz = t.z - vecEmc.Z() ;
306 0 : if(TMath::Sqrt(rawdx*rawdx+rawdz*rawdz)<fRtpc){
307 0 : Float_t dx,dz ;
308 : //calculate presize difference accounting misalignment
309 0 : GetDistanceInPHOSPlane(emcclu, t.track, dx,dz) ;
310 0 : Int_t itrack = t.track->GetID() ;
311 0 : new ((*fLinkUpArray)[iLinkUp++]) AliPHOSLink(dx, dz, iEmcRP, itrack, -1) ;
312 0 : }
313 : }
314 0 : }
315 0 : fLinkUpArray->Sort() ; //first links with smallest distances
316 0 : }
317 :
318 : //____________________________________________________________________________
319 : void AliPHOSTrackSegmentMakerv2::MakePairs()
320 : {
321 : // Using the previously made list of "links", we found the smallest link - i.e.
322 : // link with the least distance between EMC and CPV and pointing to still
323 : // unassigned RecParticles. We assign these RecPoints to TrackSegment and
324 : // remove them from the list of "unassigned".
325 :
326 : //Make arrays to mark clusters already chosen
327 : Int_t * emcExist = 0;
328 :
329 0 : if(fEmcLast > fEmcFirst)
330 0 : emcExist = new Int_t[fEmcLast-fEmcFirst] ;
331 : else
332 0 : AliFatal(Form("fEmcLast > fEmcFirst: %d > %d is not true!",fEmcLast,fEmcFirst));
333 :
334 : Int_t index;
335 0 : for(index = 0; index <fEmcLast-fEmcFirst; index ++)
336 0 : emcExist[index] = 1 ;
337 :
338 :
339 : Bool_t * tpcExist = 0;
340 0 : Int_t nTracks = fESD->GetNumberOfTracks();
341 0 : if(nTracks>0)
342 0 : tpcExist = new Bool_t[nTracks] ;
343 :
344 0 : for(index = 0; index <nTracks; index ++)
345 0 : tpcExist[index] = 1 ;
346 :
347 :
348 : // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance
349 0 : TIter nextUp(fLinkUpArray) ;
350 :
351 : AliPHOSLink * linkUp ;
352 :
353 : AliPHOSCpvRecPoint * nullpointer = 0 ;
354 :
355 0 : while ( (linkUp = static_cast<AliPHOSLink *>(nextUp()) ) ){
356 :
357 0 : if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
358 :
359 0 : if(tpcExist[linkUp->GetCpv()]){ //Track still exist
360 0 : Float_t dx,dz ;
361 0 : linkUp->GetXZ(dx,dz) ;
362 0 : new ((* fTrackSegments)[fNTrackSegments])
363 0 : AliPHOSTrackSegment(static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(linkUp->GetEmc())) ,
364 : nullpointer,
365 0 : linkUp->GetTrack(),dx,dz) ;
366 0 : (static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
367 0 : fNTrackSegments++ ;
368 0 : emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc that Cpv was found
369 : //mark track as already used
370 0 : tpcExist[linkUp->GetCpv()] = kFALSE ;
371 0 : } //if CpvUp still exist
372 : }
373 : }
374 :
375 : //look through emc recPoints left without CPV
376 0 : if(emcExist){ //if there is emc rec point
377 : Int_t iEmcRP ;
378 0 : for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst ; iEmcRP++ ){
379 0 : if(emcExist[iEmcRP] > 0 ){
380 0 : new ((*fTrackSegments)[fNTrackSegments])
381 0 : AliPHOSTrackSegment(static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP+fEmcFirst)),
382 : nullpointer) ;
383 0 : (static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
384 0 : fNTrackSegments++;
385 0 : }
386 : }
387 0 : }
388 0 : delete [] emcExist ;
389 0 : delete [] tpcExist ;
390 0 : }
391 :
392 : //____________________________________________________________________________
393 : void AliPHOSTrackSegmentMakerv2::Clusters2TrackSegments(Option_t *option)
394 : {
395 : // Steering method to perform track segment construction for events
396 : // in the range from fFirstEvent to fLastEvent.
397 : // This range is optionally set by SetEventRange().
398 : // if fLastEvent=-1 (by default), then process events until the end.
399 :
400 0 : if(strstr(option,"tim"))
401 0 : gBenchmark->Start("PHOSTSMaker");
402 :
403 0 : if(strstr(option,"print")) {
404 0 : Print() ;
405 0 : return ;
406 : }
407 :
408 : //Make some initializations
409 0 : fNTrackSegments = 0 ;
410 0 : fEmcFirst = 0 ;
411 0 : fEmcLast = 0 ;
412 :
413 0 : fTrackSegments->Clear();
414 :
415 : // if(!ReadRecPoints(ievent)) continue; //reads RecPoints for event ievent
416 :
417 0 : for(fModule = 1; fModule <= fGeom->GetNModules() ; fModule++ ) {
418 0 : FillOneModule() ;
419 0 : MakeLinks() ;
420 0 : MakePairs() ;
421 : }
422 :
423 0 : if(strstr(option,"deb"))
424 0 : PrintTrackSegments(option);
425 :
426 0 : if(strstr(option,"tim")){
427 0 : gBenchmark->Stop("PHOSTSMaker");
428 0 : Info("Exec", "took %f seconds for making TS",
429 0 : gBenchmark->GetCpuTime("PHOSTSMaker"));
430 0 : }
431 0 : }
432 :
433 : //____________________________________________________________________________
434 : void AliPHOSTrackSegmentMakerv2::Print(const Option_t *)const
435 : {
436 : // Print TrackSegmentMaker parameters
437 :
438 0 : TString message("") ;
439 0 : if( strcmp(GetName(), "") != 0 ) {
440 0 : message = "\n======== AliPHOSTrackSegmentMakerv2 ========\n" ;
441 0 : message += "Making Track segments\n" ;
442 0 : message += "with parameters:\n" ;
443 0 : message += " Maximal EMC - TPC distance (cm) %f\n" ;
444 0 : message += "============================================\n" ;
445 0 : Info("Print", message.Data(),fRtpc) ;
446 : }
447 : else
448 0 : Info("Print", "AliPHOSTrackSegmentMakerv2 not initialized ") ;
449 :
450 0 : }
451 :
452 : //____________________________________________________________________________
453 : void AliPHOSTrackSegmentMakerv2::PrintTrackSegments(Option_t * option)
454 : {
455 : // option deb - prints # of found TrackSegments
456 : // option deb all - prints as well indexed of found RecParticles assigned to the TS
457 :
458 0 : Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
459 0 : printf(" Found %d TrackSegments\n", fTrackSegments->GetEntriesFast() );
460 :
461 0 : if(strstr(option,"all")) { // printing found TS
462 0 : printf("TrackSegment # EMC RP# CPV RP#\n") ;
463 : Int_t index;
464 0 : for (index = 0 ; index <fTrackSegments->GetEntriesFast() ; index++) {
465 0 : AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )fTrackSegments->At(index) ;
466 0 : printf(" %d %d %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetTrackIndex() ) ;
467 : }
468 0 : }
469 0 : }
|
