coding convention
[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
59 ClassImp( AliPHOSTrackSegmentMakerv1) 
60
61
62 //____________________________________________________________________________
63   AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() : AliPHOSTrackSegmentMaker()
64 {
65   // default ctor (to be used mainly by Streamer)
66
67   InitParameters() ; 
68   fDefaultInit = kTRUE ; 
69 }
70
71 //____________________________________________________________________________
72  AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const TString alirunFileName, const TString eventFolderName)
73    :AliPHOSTrackSegmentMaker(alirunFileName, eventFolderName)
74 {
75   // ctor
76
77   InitParameters() ; 
78   Init() ;
79   fDefaultInit = kFALSE ; 
80 }
81
82 //____________________________________________________________________________
83  AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
84
85   // dtor
86   // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
87   if (!fDefaultInit)  
88     delete fLinkUpArray ;
89 }
90
91
92 //____________________________________________________________________________
93 const TString AliPHOSTrackSegmentMakerv1::BranchName() const 
94 {  
95  
96   return GetName() ;
97 }
98
99 //____________________________________________________________________________
100 void  AliPHOSTrackSegmentMakerv1::FillOneModule()
101 {
102   // Finds first and last indexes between which 
103   // clusters from one PHOS module are
104
105   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
106   
107   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
108   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
109  
110   //First EMC clusters
111   Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
112   for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&  
113         ((dynamic_cast<AliPHOSRecPoint *>(emcRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule ); 
114       fEmcLast ++)  ;
115   
116   //Now CPV clusters
117   Int_t totalCpv = cpvRecPoints->GetEntriesFast() ;
118
119     for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) && 
120          ((dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule ); 
121        fCpvLast ++) ;
122       
123 }
124
125 //____________________________________________________________________________
126 Float_t  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,AliPHOSRecPoint * cpvClu, Bool_t &toofar)const
127 {
128   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
129   // Clusters are sorted in "rows" and "columns" of width 1 cm
130
131   Float_t delta = 1 ;  // Width of the rows in sorting of RecPoints (in cm)
132                        // if you change this value, change it as well in xxxRecPoint::Compare()
133   Float_t r = fRcpv ;
134  
135   TVector3 vecEmc ;
136   TVector3 vecCpv ;
137   
138   emcClu->GetLocalPosition(vecEmc) ;
139   cpvClu->GetLocalPosition(vecCpv)  ; 
140
141   if(emcClu->GetPHOSMod() == cpvClu->GetPHOSMod()){ 
142     if(vecCpv.X() <= vecEmc.X() + fRcpv + 2*delta ){ 
143
144       vecCpv = vecCpv  - vecEmc ; 
145       r = vecCpv.Mag() ;
146       toofar = kFALSE ;
147
148     } // if  xPpsd >= xEmc + ...
149     else 
150       toofar = kTRUE ;
151   } 
152   else 
153     toofar = kTRUE ;
154
155   //toofar = kFALSE ;
156  
157   
158   return r ;
159 }
160
161 //____________________________________________________________________________
162 void  AliPHOSTrackSegmentMakerv1::Init()
163 {
164   // Make all memory allocations that are not possible in default constructor
165   
166   AliPHOSGetter* gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
167   
168   fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000); 
169   if ( !gime->TrackSegmentMaker() ) {
170     gime->PostTrackSegmentMaker(this);
171   }
172 }
173
174 //____________________________________________________________________________
175 void  AliPHOSTrackSegmentMakerv1::InitParameters()
176 {
177   //Initializes parameters
178   fRcpv      = 10. ;   
179   fEmcFirst  = 0 ;    
180   fEmcLast   = 0 ;   
181   fCpvFirst  = 0 ;   
182   fCpvLast   = 0 ;   
183   fLinkUpArray = 0 ;
184   fTrackSegmentsInRun       = 0 ; 
185 }
186
187
188 //____________________________________________________________________________
189 void  AliPHOSTrackSegmentMakerv1::MakeLinks()const
190
191   // Finds distances (links) between all EMC and PPSD clusters, 
192   // which are not further apart from each other than fRcpv 
193   // and sort them in accordance with this distance
194   
195   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
196   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
197   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
198
199   fLinkUpArray->Clear() ;    
200
201   AliPHOSRecPoint * cpv ;
202   AliPHOSEmcRecPoint * emcclu ;
203
204   Int_t iLinkUp  = 0 ;
205   
206   Int_t iEmcRP;
207   for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
208     emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP)) ;
209
210     Bool_t toofar ;        
211     Int_t iCpv = 0 ;    
212     for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) { 
213       
214       cpv = dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(iCpv)) ;
215       Float_t r = GetDistanceInPHOSPlane(emcclu, cpv, toofar) ;
216       
217       if(toofar)
218         break ;  
219       if(r < fRcpv) { 
220         new ((*fLinkUpArray)[iLinkUp++])  AliPHOSLink(r, iEmcRP, iCpv) ;
221       }      
222     }
223   } 
224   
225   fLinkUpArray->Sort() ;  //first links with smallest distances
226 }
227
228 //____________________________________________________________________________
229 void  AliPHOSTrackSegmentMakerv1::MakePairs()
230
231   // Using the previously made list of "links", we found the smallest link - i.e. 
232   // link with the least distance between EMC and CPV and pointing to still 
233   // unassigned RecParticles. We assign these RecPoints to TrackSegment and 
234   // remove them from the list of "unassigned". 
235
236   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
237
238   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
239   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
240   TClonesArray * trackSegments = gime->TrackSegments();
241     
242   //Make arrays to mark clusters already chosen
243   Int_t * emcExist = 0;
244   if(fEmcLast > fEmcFirst)
245     emcExist = new Int_t[fEmcLast-fEmcFirst] ;
246   
247   Int_t index;
248   for(index = 0; index <fEmcLast-fEmcFirst; index ++)
249     emcExist[index] = 1 ;
250   
251   Bool_t * cpvExist = 0;
252   if(fCpvLast > fCpvFirst)
253     cpvExist = new Bool_t[fCpvLast-fCpvFirst] ;
254   for(index = 0; index <fCpvLast-fCpvFirst; index ++)
255     cpvExist[index] = kTRUE ;
256   
257   
258   // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance 
259   TIter nextUp(fLinkUpArray) ;
260   
261   AliPHOSLink * linkUp ;
262   
263   AliPHOSRecPoint * nullpointer = 0 ;
264   
265   while ( (linkUp =  static_cast<AliPHOSLink *>(nextUp()) ) ){  
266
267     if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){ //without ppsd Up yet 
268
269       if(cpvExist[linkUp->GetPpsd()-fCpvFirst]){ //CPV still exist
270        
271        new ((* trackSegments)[fNTrackSegments]) 
272          AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(linkUp->GetEmc())) , 
273                            dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(linkUp->GetPpsd()))) ;
274        (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
275        fNTrackSegments++ ;
276        
277        emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc  that Cpv was found 
278        //mark CPV recpoint as already used 
279        cpvExist[linkUp->GetPpsd()-fCpvFirst] = kFALSE ;
280       } //if ppsdUp still exist
281     } 
282   }        
283
284   //look through emc recPoints left without CPV/PPSD
285   if(emcExist){ //if there is emc rec point
286     Int_t iEmcRP ;
287     for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst  ; iEmcRP++ ){
288       if(emcExist[iEmcRP] > 0 ){
289        new ((*trackSegments)[fNTrackSegments])  
290          AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP+fEmcFirst)), 
291                            nullpointer) ;
292        (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
293        fNTrackSegments++;    
294       } 
295     }
296   }
297   delete [] emcExist ; 
298   delete [] cpvExist ; 
299 }
300
301 //____________________________________________________________________________
302 void  AliPHOSTrackSegmentMakerv1::Exec(Option_t * option)
303 {
304   // STEERing method
305   
306   if(strstr(option,"tim"))
307     gBenchmark->Start("PHOSTSMaker");
308  
309   if(strstr(option,"print")) {
310     Print() ; 
311     return ; 
312   }
313   
314   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;  
315   
316   const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
317
318   Int_t nevents = gime->MaxEvent() ;    
319   Int_t ievent ;
320
321   for(ievent = 0; ievent < nevents; ievent++) {
322     gime->Event(ievent,"R") ;
323     //Make some initializations 
324     fNTrackSegments = 0 ;
325     fEmcFirst = 0 ;    
326     fEmcLast  = 0 ;   
327     fCpvFirst = 0 ;   
328     fCpvLast  = 0 ;   
329     
330     gime->TrackSegments()->Clear();
331
332     //    if(!ReadRecPoints(ievent))   continue; //reads RecPoints for event ievent
333     
334     for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ) {
335       FillOneModule() ; 
336       MakeLinks() ;
337       MakePairs() ;
338     }
339
340     WriteTrackSegments() ;
341
342     if(strstr(option,"deb"))
343       PrintTrackSegments(option);
344     
345     //increment the total number of track segments per run 
346     fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ; 
347
348   }
349   
350   if(strstr(option,"tim")){
351     gBenchmark->Stop("PHOSTSMaker");
352     Info("Exec", "took %f seconds for making TS %f seconds per event", 
353           gBenchmark->GetCpuTime("PHOSTSMaker"), 
354           gBenchmark->GetCpuTime("PHOSTSMaker")/nevents) ;
355    }
356   Unload();
357 }
358
359 //____________________________________________________________________________
360 void AliPHOSTrackSegmentMakerv1::Unload() 
361 {
362   // Unloads the task from the folder
363   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;  
364   gime->PhosLoader()->UnloadRecPoints() ;
365   gime->PhosLoader()->UnloadTracks() ;
366 }
367
368 //____________________________________________________________________________
369 void AliPHOSTrackSegmentMakerv1::Print()const
370 {
371   //  Print TrackSegmentMaker parameters
372
373   TString message("") ;
374   if( strcmp(GetName(), "") != 0 ) {
375     message = "\n======== AliPHOSTrackSegmentMakerv1 ========\n" ; 
376     message += "Making Track segments\n" ;
377     message += "with parameters:\n" ; 
378     message += "     Maximal EMC - CPV (PPSD) distance (cm) %f\n" ;
379     message += "============================================\n" ;
380     Info("Print", message.Data(),fRcpv) ;
381   }
382   else
383     Info("Print", "AliPHOSTrackSegmentMakerv1 not initialized ") ;
384 }
385
386 //____________________________________________________________________________
387 void AliPHOSTrackSegmentMakerv1::WriteTrackSegments()
388 {
389   // Writes found TrackSegments to TreeR. Creates branches 
390   // "PHOSTS" and "AliPHOSTrackSegmentMaker" with the same title.
391   // In the former branch found TrackSegments are stored, while 
392   // in the latter all parameters, with which TS were made. 
393   // ROOT does not allow overwriting existing branches, therefore
394   // first we check, if branches with the same title already exist.
395   // If yes - exits without writing.
396
397   AliPHOSGetter *gime = AliPHOSGetter::Instance() ; 
398
399   TClonesArray * trackSegments = gime->TrackSegments() ; 
400   trackSegments->Expand(trackSegments->GetEntriesFast()) ;
401
402   TTree * treeT = gime->TreeT();
403  
404   //First TS
405   Int_t bufferSize = 32000 ; 
406   TBranch * tsBranch = treeT->Branch("PHOSTS",&trackSegments,bufferSize);
407   tsBranch->Fill() ;  
408
409   gime->WriteTracks("OVERWRITE");
410   gime->WriteTrackSegmentMaker("OVERWRITE");
411 }
412
413
414 //____________________________________________________________________________
415 void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
416 {
417   // option deb - prints # of found TrackSegments
418   // option deb all - prints as well indexed of found RecParticles assigned to the TS
419
420   TClonesArray * trackSegments = AliPHOSGetter::Instance()->TrackSegments() ; 
421
422   Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ; 
423   printf("nevent: %d\n", gAlice->GetEvNumber()) ; 
424   printf("        Found %d TrackSegments\n", trackSegments->GetEntriesFast() ); 
425   
426   if(strstr(option,"all")) {  // printing found TS
427     printf("TrackSegment #  EMC RP#  CPV RP#\n") ; 
428     Int_t index;
429     for (index = 0 ; index <trackSegments->GetEntriesFast() ; index++) {
430       AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )trackSegments->At(index) ; 
431       printf("   %d           %d        %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetCpvIndex() ) ; 
432     }   
433   }
434 }