Change debug print-out
[u/mrichter/AliRoot.git] / PHOS / AliPHOSTrackSegmentMakerv1.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 // Implementation version 1 of algorithm class to construct PHOS track segments
18 // Track segment for PHOS is list of 
19 //        EMC RecPoint + (possibly) CPV RecPoint + (possibly) PPSD RecPoint
20 // To find TrackSegments we do the following: 
21 //  for each EMC RecPoints we look at
22 //   CPV/PPSD RecPoints in the radious fRcpv. 
23 //  If there is such a CPV RecPoint, 
24 //   we make "Link" it is just indexes of EMC and CPV/PPSD RecPoint and distance
25 //   between them in the PHOS plane. 
26 //  Then we sort "Links" and starting from the 
27 //   least "Link" pointing to the unassined EMC and CPV RecPoints assing them to 
28 //   new TrackSegment. 
29 // If there is no CPV/PPSD RecPoint we make TrackSegment 
30 // consisting from EMC alone. There is no TrackSegments without EMC RecPoint.
31 //// In principle this class should be called from AliPHOSReconstructioner, but 
32 // one can use it as well in standalone mode.
33 // Use  case:
34 //  root [0] AliPHOSTrackSegmentMakerv1 * t = new AliPHOSTrackSegmentMaker("galice.root", "tracksegmentsname", "recpointsname")
35 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
36 //               // reads gAlice from header file "galice.root", uses recpoints stored in the branch names "recpointsname" (default = "Default")
37 //               // and saves recpoints in branch named "tracksegmentsname" (default = "recpointsname")                       
38 //  root [1] t->ExecuteTask()
39 //  root [2] t->SetMaxEmcPpsdDistance(5)
40 //  root [3] t->SetTrackSegmentsBranch("max distance 5 cm")
41 //  root [4] t->ExecuteTask("deb all time") 
42 //                 
43 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH) 
44 //
45
46 // --- ROOT system ---
47 #include "TTree.h"
48 #include "TBenchmark.h"
49
50 // --- Standard library ---
51
52 // --- AliRoot header files ---
53 #include "AliPHOSGeometry.h"
54 #include "AliPHOSTrackSegmentMakerv1.h"
55 #include "AliPHOSTrackSegment.h"
56 #include "AliPHOSLink.h"
57 #include "AliPHOSGetter.h"
58 #include "AliESD.h"
59 #include "AliESDtrack.h"
60
61 ClassImp( AliPHOSTrackSegmentMakerv1) 
62
63
64 //____________________________________________________________________________
65   AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() : AliPHOSTrackSegmentMaker()
66 {
67   // default ctor (to be used mainly by Streamer)
68
69   InitParameters() ; 
70   fDefaultInit = kTRUE ; 
71 }
72
73 //____________________________________________________________________________
74  AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const TString alirunFileName, const TString eventFolderName)
75    :AliPHOSTrackSegmentMaker(alirunFileName, eventFolderName)
76 {
77   // ctor
78
79   InitParameters() ; 
80   Init() ;
81   fDefaultInit = kFALSE ; 
82   fESD = 0;
83 }
84
85 //____________________________________________________________________________
86  AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
87
88   // dtor
89   // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
90   if (!fDefaultInit)  
91     delete fLinkUpArray ;
92 }
93
94
95 //____________________________________________________________________________
96 const TString AliPHOSTrackSegmentMakerv1::BranchName() const 
97 {  
98  
99   return GetName() ;
100 }
101
102 //____________________________________________________________________________
103 void  AliPHOSTrackSegmentMakerv1::FillOneModule()
104 {
105   // Finds first and last indexes between which 
106   // clusters from one PHOS module are
107
108   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
109   
110   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
111   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
112  
113   //First EMC clusters
114   Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
115   for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&  
116         ((dynamic_cast<AliPHOSRecPoint *>(emcRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule ); 
117       fEmcLast ++)  ;
118   
119   //Now CPV clusters
120   Int_t totalCpv = cpvRecPoints->GetEntriesFast() ;
121
122     for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) && 
123          ((dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule ); 
124        fCpvLast ++) ;
125       
126 }
127
128 //____________________________________________________________________________
129 Float_t  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,AliPHOSRecPoint * cpvClu, Bool_t &toofar) const
130 {
131   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
132   // Clusters are sorted in "rows" and "columns" of width 1 cm
133
134   Float_t delta = 1 ;  // Width of the rows in sorting of RecPoints (in cm)
135                        // if you change this value, change it as well in xxxRecPoint::Compare()
136   Float_t distance = fRcpv ;
137  
138   TVector3 vecEmc ;   // Local position of EMC recpoint
139   TVector3 vecCpv ;   // Local position of CPV recpoint propagated to EMC
140   TVector3 vecDist ;  // Distance between local positions of two points
141   
142   emcClu->GetLocalPosition(vecEmc) ;
143   cpvClu->GetLocalPosition(vecCpv)  ; 
144
145   toofar = kTRUE ;
146   if(emcClu->GetPHOSMod() == cpvClu->GetPHOSMod()){ 
147
148     if (fESD != 0x0) {
149       // Extrapolate the global track direction if any to CPV and find the closest track
150       Int_t nTracks = fESD->GetNumberOfTracks();
151       Int_t iClosestTrack = -1;
152       Double_t distance;
153       Double_t minDistance = 1e6;
154       Double_t pxyz[3], xyz[3];
155       AliESDtrack *track;
156       for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
157         track = fESD->GetTrack(iTrack);
158         track->GetOuterXYZ(xyz);     // track coord on the cylinder of PHOS radius
159         track->GetOuterPxPyPz(pxyz); // track momentum ibid.
160         vecDist = PropagateToCPV(xyz,pxyz,cpvClu->GetPHOSMod());
161         vecDist -= vecCpv;
162         distance = TMath::Sqrt(vecDist.X()*vecDist.X() + vecDist.Z()*vecDist.Z());
163         // Find the closest track to the EMC recpoint
164         if (distance < minDistance) {
165           minDistance = distance;
166           iClosestTrack = iTrack;
167         }
168       }
169
170       if (iClosestTrack != -1) {
171         track = fESD->GetTrack(iClosestTrack);
172         track->GetOuterPxPyPz(pxyz); // track momentum ibid.
173         TVector3 vecCpvGlobal; // Global position of the CPV recpoint
174         AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
175         const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
176         geom->GetGlobal((AliRecPoint*)cpvClu,vecCpvGlobal);
177         for (Int_t ixyz=0; ixyz<3; ixyz++)
178           xyz[ixyz] = vecCpvGlobal[ixyz];
179         vecDist = PropagateToCPV(xyz,pxyz,cpvClu->GetPHOSMod());
180         vecDist -= vecEmc;
181         distance = TMath::Sqrt(vecDist.X()*vecDist.X() + vecDist.Z()*vecDist.Z());
182       }
183     } else {
184       // If no ESD, than simply find EMC-CPV distance
185       distance = (vecCpv - vecEmc).Mag() ;
186     }
187     if(distance < fRcpv + 2*delta )
188       toofar = kFALSE ;
189   } 
190   
191   return distance ;
192 }
193
194 //____________________________________________________________________________
195 TVector3  AliPHOSTrackSegmentMakerv1::PropagateToCPV(Double_t *x, Double_t *p,
196                                                      Int_t moduleNumber) const
197 {
198   // Propagate a straight-line track from the origin point x
199   // along the direction p to the CPV module moduleNumber
200   // Returns a local position of such a propagation
201
202   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
203   const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
204   TVector3 moduleCenter = geom->GetCpvModuleCenter(moduleNumber);
205   TVector3 vertex(x);
206   TVector3 direction(p);
207
208   Double_t time = (moduleCenter.Mag2() - vertex.Dot(moduleCenter)) /
209     (direction.Dot(moduleCenter));
210   TVector3 globalIntersection = vertex + direction*time;
211   return geom->Global2LocalCpv(globalIntersection,moduleNumber);
212 }
213
214 //____________________________________________________________________________
215 void  AliPHOSTrackSegmentMakerv1::Init()
216 {
217   // Make all memory allocations that are not possible in default constructor
218   
219   AliPHOSGetter* gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
220   
221   fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000); 
222   if ( !gime->TrackSegmentMaker() ) {
223     gime->PostTrackSegmentMaker(this);
224   }
225 }
226
227 //____________________________________________________________________________
228 void  AliPHOSTrackSegmentMakerv1::InitParameters()
229 {
230   //Initializes parameters
231   fRcpv      = 10. ;   
232   fEmcFirst  = 0 ;    
233   fEmcLast   = 0 ;   
234   fCpvFirst  = 0 ;   
235   fCpvLast   = 0 ;   
236   fLinkUpArray = 0 ;
237   fTrackSegmentsInRun       = 0 ; 
238   SetEventRange(0,-1) ;
239 }
240
241
242 //____________________________________________________________________________
243 void  AliPHOSTrackSegmentMakerv1::MakeLinks()const
244
245   // Finds distances (links) between all EMC and PPSD clusters, 
246   // which are not further apart from each other than fRcpv 
247   // and sort them in accordance with this distance
248   
249   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
250   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
251   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
252
253   fLinkUpArray->Clear() ;    
254
255   AliPHOSRecPoint * cpv ;
256   AliPHOSEmcRecPoint * emcclu ;
257
258   Int_t iLinkUp  = 0 ;
259   
260   Int_t iEmcRP;
261   for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
262     emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP)) ;
263
264     Bool_t toofar ;        
265     Int_t iCpv = 0 ;    
266     for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) { 
267       
268       cpv = dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(iCpv)) ;
269       Float_t r = GetDistanceInPHOSPlane(emcclu, cpv, toofar) ;
270       
271       if(toofar)
272         break ;  
273       if(r < fRcpv) { 
274         new ((*fLinkUpArray)[iLinkUp++])  AliPHOSLink(r, iEmcRP, iCpv) ;
275       }      
276     }
277   } 
278   
279   fLinkUpArray->Sort() ;  //first links with smallest distances
280 }
281
282 //____________________________________________________________________________
283 void  AliPHOSTrackSegmentMakerv1::MakePairs()
284
285   // Using the previously made list of "links", we found the smallest link - i.e. 
286   // link with the least distance between EMC and CPV and pointing to still 
287   // unassigned RecParticles. We assign these RecPoints to TrackSegment and 
288   // remove them from the list of "unassigned". 
289
290   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
291
292   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
293   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
294   TClonesArray * trackSegments = gime->TrackSegments();
295     
296   //Make arrays to mark clusters already chosen
297   Int_t * emcExist = 0;
298   if(fEmcLast > fEmcFirst)
299     emcExist = new Int_t[fEmcLast-fEmcFirst] ;
300   
301   Int_t index;
302   for(index = 0; index <fEmcLast-fEmcFirst; index ++)
303     emcExist[index] = 1 ;
304   
305   Bool_t * cpvExist = 0;
306   if(fCpvLast > fCpvFirst)
307     cpvExist = new Bool_t[fCpvLast-fCpvFirst] ;
308   for(index = 0; index <fCpvLast-fCpvFirst; index ++)
309     cpvExist[index] = kTRUE ;
310   
311   
312   // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance 
313   TIter nextUp(fLinkUpArray) ;
314   
315   AliPHOSLink * linkUp ;
316   
317   AliPHOSRecPoint * nullpointer = 0 ;
318   
319   while ( (linkUp =  static_cast<AliPHOSLink *>(nextUp()) ) ){  
320
321     if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){ //without ppsd Up yet 
322
323       if(cpvExist[linkUp->GetPpsd()-fCpvFirst]){ //CPV still exist
324        
325        new ((* trackSegments)[fNTrackSegments]) 
326          AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(linkUp->GetEmc())) , 
327                            dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(linkUp->GetPpsd()))) ;
328        (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
329        fNTrackSegments++ ;
330        
331        emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc  that Cpv was found 
332        //mark CPV recpoint as already used 
333        cpvExist[linkUp->GetPpsd()-fCpvFirst] = kFALSE ;
334       } //if ppsdUp still exist
335     } 
336   }        
337
338   //look through emc recPoints left without CPV/PPSD
339   if(emcExist){ //if there is emc rec point
340     Int_t iEmcRP ;
341     for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst  ; iEmcRP++ ){
342       if(emcExist[iEmcRP] > 0 ){
343        new ((*trackSegments)[fNTrackSegments])  
344          AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP+fEmcFirst)), 
345                            nullpointer) ;
346        (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
347        fNTrackSegments++;    
348       } 
349     }
350   }
351   delete [] emcExist ; 
352   delete [] cpvExist ; 
353 }
354
355 //____________________________________________________________________________
356 void  AliPHOSTrackSegmentMakerv1::Exec(Option_t *option)
357 {
358   // Steering method to perform track segment construction for events
359   // in the range from fFirstEvent to fLastEvent.
360   // This range is optionally set by SetEventRange().
361   // if fLastEvent=-1 (by default), then process events until the end.
362   
363   if(strstr(option,"tim"))
364     gBenchmark->Start("PHOSTSMaker");
365  
366   if(strstr(option,"print")) {
367     Print() ; 
368     return ; 
369   }
370   
371   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ;  
372  
373   const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
374
375   if (fLastEvent == -1) fLastEvent = gime->MaxEvent() - 1 ;
376   else fLastEvent = TMath::Min(fFirstEvent,gime->MaxEvent());
377   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
378
379   Int_t ievent ; 
380   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
381     gime->Event(ievent,"R") ;
382     //Make some initializations 
383     fNTrackSegments = 0 ;
384     fEmcFirst = 0 ;    
385     fEmcLast  = 0 ;   
386     fCpvFirst = 0 ;   
387     fCpvLast  = 0 ;   
388     
389     gime->TrackSegments()->Clear();
390
391     //    if(!ReadRecPoints(ievent))   continue; //reads RecPoints for event ievent
392     
393     for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ) {
394       FillOneModule() ; 
395       MakeLinks() ;
396       MakePairs() ;
397     }
398
399     WriteTrackSegments() ;
400
401     if(strstr(option,"deb"))
402       PrintTrackSegments(option);
403     
404     //increment the total number of track segments per run 
405     fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ; 
406
407   }
408   
409   if(strstr(option,"tim")){
410     gBenchmark->Stop("PHOSTSMaker");
411     Info("Exec", "took %f seconds for making TS %f seconds per event", 
412           gBenchmark->GetCpuTime("PHOSTSMaker"), 
413           gBenchmark->GetCpuTime("PHOSTSMaker")/nEvents) ;
414    }
415   Unload();
416 }
417
418 //____________________________________________________________________________
419 void AliPHOSTrackSegmentMakerv1::Unload() 
420 {
421   // Unloads the task from the folder
422   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;  
423   gime->PhosLoader()->UnloadRecPoints() ;
424   gime->PhosLoader()->UnloadTracks() ;
425 }
426
427 //____________________________________________________________________________
428 void AliPHOSTrackSegmentMakerv1::Print()const
429 {
430   //  Print TrackSegmentMaker parameters
431
432   TString message("") ;
433   if( strcmp(GetName(), "") != 0 ) {
434     message = "\n======== AliPHOSTrackSegmentMakerv1 ========\n" ; 
435     message += "Making Track segments\n" ;
436     message += "with parameters:\n" ; 
437     message += "     Maximal EMC - CPV (PPSD) distance (cm) %f\n" ;
438     message += "============================================\n" ;
439     Info("Print", message.Data(),fRcpv) ;
440   }
441   else
442     Info("Print", "AliPHOSTrackSegmentMakerv1 not initialized ") ;
443 }
444
445 //____________________________________________________________________________
446 void AliPHOSTrackSegmentMakerv1::WriteTrackSegments()
447 {
448   // Writes found TrackSegments to TreeR. Creates branches 
449   // "PHOSTS" and "AliPHOSTrackSegmentMaker" with the same title.
450   // In the former branch found TrackSegments are stored, while 
451   // in the latter all parameters, with which TS were made. 
452   // ROOT does not allow overwriting existing branches, therefore
453   // first we check, if branches with the same title already exist.
454   // If yes - exits without writing.
455
456   AliPHOSGetter *gime = AliPHOSGetter::Instance() ; 
457
458   TClonesArray * trackSegments = gime->TrackSegments() ; 
459   trackSegments->Expand(trackSegments->GetEntriesFast()) ;
460
461   TTree * treeT = gime->TreeT();
462  
463   //First TS
464   Int_t bufferSize = 32000 ; 
465   TBranch * tsBranch = treeT->Branch("PHOSTS",&trackSegments,bufferSize);
466   tsBranch->Fill() ;  
467
468   gime->WriteTracks("OVERWRITE");
469   gime->WriteTrackSegmentMaker("OVERWRITE");
470 }
471
472
473 //____________________________________________________________________________
474 void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
475 {
476   // option deb - prints # of found TrackSegments
477   // option deb all - prints as well indexed of found RecParticles assigned to the TS
478
479   TClonesArray * trackSegments = AliPHOSGetter::Instance()->TrackSegments() ; 
480
481   Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ; 
482   printf("nevent: %d\n", gAlice->GetEvNumber()) ; 
483   printf("        Found %d TrackSegments\n", trackSegments->GetEntriesFast() ); 
484   
485   if(strstr(option,"all")) {  // printing found TS
486     printf("TrackSegment #  EMC RP#  CPV RP#\n") ; 
487     Int_t index;
488     for (index = 0 ; index <trackSegments->GetEntriesFast() ; index++) {
489       AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )trackSegments->At(index) ; 
490       printf("   %d           %d        %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetCpvIndex() ) ; 
491     }   
492   }
493 }