Added README to PHOS/macros/Trigger/OCDB, explaining the content.
[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)
6ba1dd81 104
105 for(Int_t i=0; i<5; i++)fNtpcTracks[i]=0 ;
861b23c3 106 InitParameters() ;
861b23c3 107}
108
109//____________________________________________________________________________
9a2cdbdf 110AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(AliPHOSGeometry *geom) :
111 AliPHOSTrackSegmentMaker(geom),
861b23c3 112 fDefaultInit(kFALSE),
113 fWrite(kFALSE),
114 fNTrackSegments(0),
115 fRtpc(0.f),
c978a69b 116 fVtx(0.,0.,0.),
861b23c3 117 fLinkUpArray(0),
c978a69b 118 fTPCtracks(),
861b23c3 119 fEmcFirst(0),
120 fEmcLast(0),
121 fModule(0),
9a2cdbdf 122 fTrackSegments(NULL)
861b23c3 123{
124 // ctor
6ba1dd81 125 for(Int_t i=0; i<5; i++)fNtpcTracks[i]=0 ;
126
861b23c3 127 InitParameters() ;
128 Init() ;
129 fESD = 0;
130}
131
132//____________________________________________________________________________
133AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const AliPHOSTrackSegmentMakerv2 & tsm) :
134 AliPHOSTrackSegmentMaker(tsm),
135 fDefaultInit(kFALSE),
136 fWrite(kFALSE),
137 fNTrackSegments(0),
138 fRtpc(0.f),
c978a69b 139 fVtx(0.,0.,0.),
861b23c3 140 fLinkUpArray(0),
c978a69b 141 fTPCtracks(),
861b23c3 142 fEmcFirst(0),
143 fEmcLast(0),
144 fModule(0),
9a2cdbdf 145 fTrackSegments(NULL)
861b23c3 146{
147 // cpy ctor: no implementation yet
148 // requested by the Coding Convention
149 Fatal("cpy ctor", "not implemented") ;
150}
151
152
153//____________________________________________________________________________
154 AliPHOSTrackSegmentMakerv2::~AliPHOSTrackSegmentMakerv2()
155{
156 // dtor
157 // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
158 if (!fDefaultInit)
159 delete fLinkUpArray ;
9a2cdbdf 160 if (fTrackSegments) {
161 fTrackSegments->Delete();
162 delete fTrackSegments;
163 }
861b23c3 164}
165
166//____________________________________________________________________________
167void AliPHOSTrackSegmentMakerv2::FillOneModule()
168{
169 // Finds first and last indexes between which
170 // clusters from one PHOS module are
171
861b23c3 172 //First EMC clusters
9a2cdbdf 173 Int_t totalEmc = fEMCRecPoints->GetEntriesFast() ;
861b23c3 174 for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
c78764fb 175 ((static_cast<AliPHOSRecPoint *>(fEMCRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule );
861b23c3 176 fEmcLast ++) ;
177
178 //Now TPC tracks
4b598bb1 179 if(fESD){
180 //Do it ones, only first time
181 if(fModule==1){
182 Int_t nTracks = fESD->GetNumberOfTracks();
c978a69b 183
9a2cdbdf 184 Int_t nPHOSmod = fGeom->GetNModules() ;
c978a69b 185 if(fTPCtracks[0].size()<(UInt_t)nTracks){
186 for(Int_t i=0; i<nPHOSmod; i++)
187 fTPCtracks[i].resize(nTracks) ;
188 }
189 for(Int_t i=0; i<5; i++)fNtpcTracks[i]=0 ;
4b598bb1 190 TVector3 inPHOS ;
191
192 //In this particular case we use fixed vertex position at zero
193 Double_t vtx[3]={0.,0.,0.} ;
194 AliESDtrack *track;
195 Double_t xyz[3] ;
9a2cdbdf 196 Double_t rEMC = fGeom->GetIPtoCrystalSurface() ; //Use here ideal geometry
4b598bb1 197 for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
198 track = fESD->GetTrack(iTrack);
c978a69b 199 for(Int_t iTestMod=1; iTestMod<=nPHOSmod; iTestMod++){
9a2cdbdf 200 Double_t modAngle=270.+fGeom->GetPHOSAngle(iTestMod) ;
c978a69b 201 modAngle*=TMath::Pi()/180. ;
202 track->Rotate(modAngle);
203 if (!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz))
204 continue; //track coord on the cylinder of PHOS radius
205 if ((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0)
206 continue;
207 //Check if this track hits PHOS
208 inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
209 Int_t modNum ;
210 Double_t x,z ;
9a2cdbdf 211 fGeom->ImpactOnEmc(vtx, inPHOS.Theta(), inPHOS.Phi(), modNum, z, x) ;
c978a69b 212 if(modNum==iTestMod){
213 //Mark this track as one belonging to module
214 TrackInPHOS_t &t = fTPCtracks[modNum-1][fNtpcTracks[modNum-1]++] ;
215 t.track = track ;
216 t.x = x ;
217 t.z = z ;
218 break ;
219 }
4b598bb1 220 }
221 }
861b23c3 222 }
223 }
4b598bb1 224
861b23c3 225}
226
227//____________________________________________________________________________
228void 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
c978a69b 235 TVector3 emcGlobal; // Global position of the EMC recpoint
9a2cdbdf 236 fGeom->GetGlobalPHOS((AliPHOSRecPoint*)emcClu,emcGlobal);
4b598bb1 237
c978a69b 238 //Calculate actual distance to the PHOS surface
9a2cdbdf 239 Double_t modAngle=270.+fGeom->GetPHOSAngle(emcClu->GetPHOSMod()) ;
c978a69b 240 modAngle*=TMath::Pi()/180. ;
241 Double_t rEMC = emcGlobal.X()*TMath::Cos(modAngle)+emcGlobal.Y()*TMath::Sin(modAngle) ;
242 track->Rotate(modAngle);
861b23c3 243 Double_t xyz[3] ;
c978a69b 244 if(!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz)){
245 dx=999. ;
246 dz=999. ;
247 return ; //track coord on the cylinder of PHOS radius
861b23c3 248 }
c978a69b 249 if((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0){
861b23c3 250 dx=999. ;
251 dz=999. ;
c978a69b 252 return ;
253 }
254
255 dx=(emcGlobal.X()-xyz[0])*TMath::Sin(modAngle)-(emcGlobal.Y()-xyz[1])*TMath::Cos(modAngle) ;
256 dx=TMath::Sign(dx,(Float_t)(emcGlobal.X()-xyz[0])) ; //set direction
257 dz=emcGlobal.Z()-xyz[2] ;
861b23c3 258
259 return ;
260}
261//____________________________________________________________________________
262void AliPHOSTrackSegmentMakerv2::Init()
263{
264 // Make all memory allocations that are not possible in default constructor
265
861b23c3 266 fLinkUpArray = new TClonesArray("AliPHOSLink", 1000);
9a2cdbdf 267 fTrackSegments = new TClonesArray("AliPHOSTrackSegment",100);
861b23c3 268}
269
270//____________________________________________________________________________
271void AliPHOSTrackSegmentMakerv2::InitParameters()
272{
273 //Initializes parameters
274 fRtpc = 4. ;
275 fEmcFirst = 0 ;
276 fEmcLast = 0 ;
277 fLinkUpArray = 0 ;
278 fWrite = kTRUE ;
861b23c3 279}
280
281
282//____________________________________________________________________________
c978a69b 283void AliPHOSTrackSegmentMakerv2::MakeLinks()
861b23c3 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
861b23c3 289 fLinkUpArray->Clear() ;
290
291 AliPHOSEmcRecPoint * emcclu ;
292
293 Int_t iLinkUp = 0 ;
294
295 Int_t iEmcRP;
296 for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
c78764fb 297 emcclu = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP)) ;
c978a69b 298 TVector3 vecEmc ;
299 emcclu->GetLocalPosition(vecEmc) ;
861b23c3 300 Int_t mod=emcclu->GetPHOSMod() ;
c978a69b 301 for(Int_t itr=0; itr<fNtpcTracks[mod-1] ; itr++){
302 TrackInPHOS_t &t = fTPCtracks[mod-1][itr] ;
303 //calculate raw distance
304 Double_t rawdx = t.x - vecEmc.X() ;
305 Double_t rawdz = t.z - vecEmc.Z() ;
306 if(TMath::Sqrt(rawdx*rawdx+rawdz*rawdz)<fRtpc){
307 Float_t dx,dz ;
308 //calculate presize difference accounting misalignment
309 GetDistanceInPHOSPlane(emcclu, t.track, dx,dz) ;
310 Int_t itrack = t.track->GetID() ;
861b23c3 311 new ((*fLinkUpArray)[iLinkUp++]) AliPHOSLink(dx, dz, iEmcRP, itrack, -1) ;
312 }
313 }
314 }
861b23c3 315 fLinkUpArray->Sort() ; //first links with smallest distances
316}
317
318//____________________________________________________________________________
319void 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
861b23c3 326 //Make arrays to mark clusters already chosen
327 Int_t * emcExist = 0;
c78764fb 328
861b23c3 329 if(fEmcLast > fEmcFirst)
330 emcExist = new Int_t[fEmcLast-fEmcFirst] ;
c78764fb 331 else
332 AliFatal(Form("fEmcLast > fEmcFirst: %d > %d is not true!",fEmcLast,fEmcFirst));
861b23c3 333
334 Int_t index;
335 for(index = 0; index <fEmcLast-fEmcFirst; index ++)
336 emcExist[index] = 1 ;
337
338
c978a69b 339 Bool_t * tpcExist = 0;
340 Int_t nTracks = fESD->GetNumberOfTracks();
861b23c3 341 if(nTracks>0)
c978a69b 342 tpcExist = new Bool_t[nTracks] ;
861b23c3 343
344 for(index = 0; index <nTracks; index ++)
345 tpcExist[index] = 1 ;
346
347
348 // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance
349 TIter nextUp(fLinkUpArray) ;
350
351 AliPHOSLink * linkUp ;
352
353 AliPHOSCpvRecPoint * nullpointer = 0 ;
354
355 while ( (linkUp = static_cast<AliPHOSLink *>(nextUp()) ) ){
356
357 if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
358
359 if(tpcExist[linkUp->GetCpv()]){ //Track still exist
360 Float_t dx,dz ;
361 linkUp->GetXZ(dx,dz) ;
9a2cdbdf 362 new ((* fTrackSegments)[fNTrackSegments])
c78764fb 363 AliPHOSTrackSegment(static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(linkUp->GetEmc())) ,
861b23c3 364 nullpointer,
365 linkUp->GetTrack(),dx,dz) ;
c78764fb 366 (static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
861b23c3 367 fNTrackSegments++ ;
368 emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc that Cpv was found
369 //mark track as already used
370 tpcExist[linkUp->GetCpv()] = kFALSE ;
371 } //if CpvUp still exist
372 }
373 }
374
375 //look through emc recPoints left without CPV
376 if(emcExist){ //if there is emc rec point
377 Int_t iEmcRP ;
378 for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst ; iEmcRP++ ){
379 if(emcExist[iEmcRP] > 0 ){
9a2cdbdf 380 new ((*fTrackSegments)[fNTrackSegments])
c78764fb 381 AliPHOSTrackSegment(static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP+fEmcFirst)),
861b23c3 382 nullpointer) ;
97e75f77 383 (static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
861b23c3 384 fNTrackSegments++;
385 }
386 }
387 }
388 delete [] emcExist ;
389 delete [] tpcExist ;
390}
391
392//____________________________________________________________________________
9a2cdbdf 393void AliPHOSTrackSegmentMakerv2::Clusters2TrackSegments(Option_t *option)
861b23c3 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 if(strstr(option,"tim"))
401 gBenchmark->Start("PHOSTSMaker");
402
403 if(strstr(option,"print")) {
404 Print() ;
405 return ;
406 }
407
9a2cdbdf 408 //Make some initializations
409 fNTrackSegments = 0 ;
410 fEmcFirst = 0 ;
411 fEmcLast = 0 ;
861b23c3 412
9a2cdbdf 413 fTrackSegments->Clear();
861b23c3 414
9a2cdbdf 415 // if(!ReadRecPoints(ievent)) continue; //reads RecPoints for event ievent
861b23c3 416
9a2cdbdf 417 for(fModule = 1; fModule <= fGeom->GetNModules() ; fModule++ ) {
418 FillOneModule() ;
419 MakeLinks() ;
420 MakePairs() ;
421 }
861b23c3 422
9a2cdbdf 423 if(strstr(option,"deb"))
424 PrintTrackSegments(option);
861b23c3 425
861b23c3 426 if(strstr(option,"tim")){
427 gBenchmark->Stop("PHOSTSMaker");
9a2cdbdf 428 Info("Exec", "took %f seconds for making TS",
429 gBenchmark->GetCpuTime("PHOSTSMaker"));
430 }
861b23c3 431}
432
433//____________________________________________________________________________
434void AliPHOSTrackSegmentMakerv2::Print(const Option_t *)const
435{
436 // Print TrackSegmentMaker parameters
437
438 TString message("") ;
439 if( strcmp(GetName(), "") != 0 ) {
440 message = "\n======== AliPHOSTrackSegmentMakerv2 ========\n" ;
441 message += "Making Track segments\n" ;
442 message += "with parameters:\n" ;
443 message += " Maximal EMC - TPC distance (cm) %f\n" ;
444 message += "============================================\n" ;
445 Info("Print", message.Data(),fRtpc) ;
446 }
447 else
448 Info("Print", "AliPHOSTrackSegmentMakerv2 not initialized ") ;
861b23c3 449
861b23c3 450}
451
861b23c3 452//____________________________________________________________________________
453void 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
861b23c3 458 Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
9a2cdbdf 459 printf(" Found %d TrackSegments\n", fTrackSegments->GetEntriesFast() );
861b23c3 460
461 if(strstr(option,"all")) { // printing found TS
462 printf("TrackSegment # EMC RP# CPV RP#\n") ;
463 Int_t index;
9a2cdbdf 464 for (index = 0 ; index <fTrackSegments->GetEntriesFast() ; index++) {
465 AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )fTrackSegments->At(index) ;
861b23c3 466 printf(" %d %d %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetTrackIndex() ) ;
467 }
468 }
469}