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