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