Transition to NewIO
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALTrackSegmentMakerv1.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 EMCAL track segments
18 // Track segment for EMCAL is list of 
19 //        ECAL RecPoint + (possibly) PRE RecPoint + (possibly) HCAL RecPoint
20 // To find TrackSegments we do the following: 
21 //  for each ECAL RecPoint we look for PRE and HC RecPoints with same direction  within fSame. 
22 //  If there is such a PRE or ECAL RecPoint, 
23 //   we make a "Link": indexes of ECAL and PRE, HCAL  RecPoints and their scalar product. 
24 //  Then we sort "Links", starting from the 
25 //   least "Link" pointing to the unassigned RecPoints assigning them to a new TrackSegment. 
26 //  If there is no PRE, HCAL RecPoint we make a TrackSegment 
27 //   consisting from ECAL alone. There is no TrackSegments without ECAL RecPoint.
28 //// In principle this class should be called from AliEMCALReconstructioner, but 
29 // one can use it as well in standalone mode.
30 //                 
31 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH) 
32 //
33
34 // --- ROOT system ---
35 #include "TROOT.h"
36 #include "TFile.h"
37 #include "TFolder.h"
38 #include "TTree.h"
39 #include "TSystem.h"
40 #include "TBenchmark.h"
41
42 // --- Standard library ---
43
44 // --- AliRoot header files ---
45
46 #include "AliEMCALTrackSegmentMakerv1.h"
47 #include "AliEMCALClusterizerv1.h"
48 #include "AliEMCALTrackSegment.h"
49 #include "AliEMCALLink.h"
50 #include "AliEMCALGetter.h"
51 #include "AliEMCAL.h"
52
53 ClassImp( AliEMCALTrackSegmentMakerv1) 
54
55
56 //____________________________________________________________________________
57   AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1() : AliEMCALTrackSegmentMaker()
58 {
59   // default ctor (to be used mainly by Streamer)
60
61   InitParameters() ; 
62   fDefaultInit = kTRUE ; 
63 }
64
65 //____________________________________________________________________________
66  AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1(const TString alirunFileName, const TString eventFolderName)
67    :AliEMCALTrackSegmentMaker(alirunFileName, eventFolderName)
68 {
69   // ctor
70
71   InitParameters() ; 
72   Init() ;
73   fDefaultInit = kFALSE ; 
74
75 }
76
77 //____________________________________________________________________________
78  AliEMCALTrackSegmentMakerv1::~AliEMCALTrackSegmentMakerv1()
79
80   // dtor
81   // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
82   
83   if (!fDefaultInit) {
84     delete fPRELinkArray  ;
85     delete fHCALinkArray  ;
86   }
87 }
88
89 //____________________________________________________________________________
90 const TString AliEMCALTrackSegmentMakerv1::BranchName() const 
91 {  
92    return GetName() ;
93
94 }
95
96 //____________________________________________________________________________
97 Float_t  AliEMCALTrackSegmentMakerv1::HowClose(AliEMCALTowerRecPoint * ec, AliEMCALTowerRecPoint * rp, Bool_t &toofar)const
98 {
99   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
100   // Clusters are sorted in "rows" and "columns" of width 1 cm
101
102   Float_t r = -1. ;
103   Float_t delta = 10. ;  // large enough to be ineffective ??! 
104
105  
106   TVector3 vecEC = ec->XYZInAlice() ;
107   TVector3 vecRP = rp->XYZInAlice() ;
108   
109   Float_t pro = TMath::Abs(1 - (vecEC * vecRP / ( vecEC.Mag() * vecRP.Mag() ))) ; 
110
111   if(pro <= delta) { 
112     r = pro ;
113     toofar = kFALSE ;
114   } 
115   else 
116     toofar = kTRUE ;
117
118   if (gDebug == 2 ) 
119     Info("HowClose", "ec = %d, rp = %d pro = %f, toofar=%d", ec->GetIndexInList(), rp->GetIndexInList(), pro, toofar ) ; 
120  
121   return r ;
122 }
123
124 //____________________________________________________________________________
125 void  AliEMCALTrackSegmentMakerv1::Init()
126 {
127   // Make all memory allocations that are not possible in default constructor
128   
129   AliEMCALGetter* gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());
130
131   fPRELinkArray = new TClonesArray("AliEMCALLink", 1000); 
132   fHCALinkArray = new TClonesArray("AliEMCALLink", 1000); 
133   if ( !gime->TrackSegmentMaker() ) {
134     gime->PostTrackSegmentMaker(this);
135   }
136 }
137
138 //____________________________________________________________________________
139 void  AliEMCALTrackSegmentMakerv1::InitParameters()
140 {
141   fClose              = 10e-3 ;   
142   fPRELinkArray       = 0 ;
143   fHCALinkArray       = 0 ;
144   fTrackSegmentsInRun = 0 ; 
145 }
146
147
148 //____________________________________________________________________________
149 void  AliEMCALTrackSegmentMakerv1::MakeLinks()const
150
151   // Finds distances (links) between all PRE, EC and HC clusters, 
152   // which are not further apart from each other than fDangle 
153   // and sort them in accordance with this distance
154   
155   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
156   TObjArray * aECARecPoints  = gime->ECARecPoints() ; 
157   TObjArray * aPRERecPoints  = gime->PRERecPoints() ; 
158   TObjArray * aHCARecPoints  = gime->HCARecPoints() ; 
159
160   fPRELinkArray->Clear() ;    
161   fHCALinkArray->Clear() ;    
162
163   AliEMCALTowerRecPoint * pre ;
164   AliEMCALTowerRecPoint * eca ;
165   AliEMCALTowerRecPoint * hca ;
166
167   Int_t iPRELink  = 0 ;
168   Int_t iHCALink  = 0 ;
169     
170   Int_t iECARP;
171   for(iECARP = 0; iECARP < aECARecPoints->GetEntriesFast(); iECARP++ ) {
172     eca = dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(iECARP)) ;
173     Bool_t toofar = kTRUE ;        
174     Int_t iPRERP = 0 ;    
175     for(iPRERP = 0; iPRERP < aPRERecPoints->GetEntriesFast(); iPRERP++ ) { 
176       pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(iPRERP)) ;
177       Float_t prod = HowClose(eca, pre, toofar) ;    
178       if(toofar)
179         break ;  
180       if(prod < fClose) { 
181         new ((*fPRELinkArray)[iPRELink++])  AliEMCALLink(prod, iECARP, iPRERP, 0) ;
182       }      
183     }
184     toofar = kTRUE ; 
185     Int_t iHCARP = 0 ;    
186     for(iHCARP = 0; iHCARP < aHCARecPoints->GetEntriesFast(); iHCARP++ ) { 
187       hca = dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(iHCARP)) ;
188       Float_t prod = HowClose(eca, hca, toofar) ;    
189       if(toofar)
190         break ;  
191       if(prod < fClose) { 
192         new ((*fHCALinkArray)[iHCALink++])  AliEMCALLink(prod, iECARP, iHCARP, 1) ;
193       }      
194     }
195   } 
196   
197   fPRELinkArray->Sort() ;  //first links with largest scalar product
198   fHCALinkArray->Sort() ;  //first links with largest scalar product
199 }
200
201 //____________________________________________________________________________
202 void  AliEMCALTrackSegmentMakerv1::MakePairs()
203
204   // Using the previously made list of "links", we found the best link - i.e. 
205   // link with the largest scalar product (closest to one) to still 
206   // unassigned RecParticles. We assign these RecPoints to TrackSegment and 
207   // remove them from the list of "unassigned". 
208   
209   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
210   TObjArray * aECARecPoints = gime->ECARecPoints() ; 
211   TObjArray * aPRERecPoints = gime->PRERecPoints() ; 
212   TObjArray * aHCARecPoints = gime->HCARecPoints() ; 
213   TClonesArray * trackSegments = gime->TrackSegments() ;   
214     
215   //Make arrays to mark clusters already chosen
216   Int_t * ecaExist = 0;
217   Int_t nECA = aECARecPoints->GetEntriesFast() ;  
218   if (nECA) 
219     ecaExist = new Int_t[nECA] ;
220   
221   Int_t index;
222   for(index = 0; index < nECA; index ++)
223     ecaExist[index] = 1 ;
224   
225   Bool_t * preExist = 0;
226   Int_t nPRE = aPRERecPoints->GetEntriesFast() ;  
227   if(nPRE)
228     preExist = new Bool_t[nPRE] ;
229   for(index = 0; index < nPRE; index ++)
230     preExist[index] = kTRUE ;
231   
232   Bool_t * hcaExist = 0;
233   Int_t nHCA = aHCARecPoints->GetEntriesFast() ;  
234   if(nHCA)
235     hcaExist = new Bool_t[nHCA] ;
236   for(index = 0; index < nHCA; index ++)
237     hcaExist[index] = kTRUE ;
238   
239   AliEMCALTowerRecPoint * null = 0 ; 
240  // Finds the smallest links and makes pairs of PRE and ECAL clusters with largest scalar product 
241  
242   TIter nextPRE(fPRELinkArray) ;
243   AliEMCALLink * linkPRE ;
244   
245   while ( (linkPRE =  static_cast<AliEMCALLink *>(nextPRE()) ) ){  
246
247     if(ecaExist[linkPRE->GetECA()] != -1){ //without PRE yet 
248
249       if(preExist[linkPRE->GetOther()]){ // PRE still exist
250         
251         new ((* trackSegments)[fNTrackSegments]) 
252           AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(linkPRE->GetECA())) , 
253                                dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(linkPRE->GetOther())), null) ;
254         (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
255         fNTrackSegments++ ;
256         if (gDebug == 2 ) 
257           Info("MakePairs", "ECAL section with PRE section") ;  
258         ecaExist[linkPRE->GetECA()] = -1 ; //Mark ecal  that pre was found 
259         //mark PRE recpoint as already used 
260         preExist[linkPRE->GetOther()] = kFALSE ;
261       } //if PRE still exist
262     } 
263   }
264
265   // Finds the smallest links and makes pairs of HCAL and ECAL clusters with largest scalar product 
266  
267   TIter nextHCA(fHCALinkArray) ;
268   AliEMCALLink * linkHCA ;
269   
270   while ( (linkHCA =  static_cast<AliEMCALLink *>(nextHCA()) ) ){  
271
272     if(ecaExist[linkHCA->GetECA()] != -2){ //without HCAL yet 
273
274       if(hcaExist[linkHCA->GetOther()]){ // HCAL still exist
275         // search among the already existing track segments
276         Int_t ii ; 
277         Bool_t found = kFALSE ; 
278         AliEMCALTrackSegment * ts = 0 ; 
279         for ( ii = 0 ; ii < fNTrackSegments ; ii++ ) {
280           ts = dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(ii)) ;
281           if ( ts->GetECAIndex() == linkHCA->GetECA() ) {
282             found = kTRUE ; 
283             break ; 
284           }
285         }
286         if (found){
287           ts->SetHCARecPoint( dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(linkHCA->GetOther())) ) ;
288           if (gDebug == 2 ) 
289             Info("MakePairs", "ECAL section with PRE and HCAL sections") ;
290         }       
291         if (!found) {
292           new ((* trackSegments)[fNTrackSegments]) 
293             AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(linkHCA->GetECA())), null,
294                                  dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(linkHCA->GetOther()))) ; 
295         (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
296         fNTrackSegments++ ;
297         if (gDebug == 2 ) 
298           Info("MakePairs", "ECAL section with HCAL section") ;         
299         }
300         ecaExist[linkHCA->GetECA()] = -2 ; //Mark ecal  that hcal was found 
301         //mark HCAL recpoint as already used 
302         hcaExist[linkHCA->GetOther()] = kFALSE ;
303       } //if HCAL still exist
304     } 
305   }
306   
307
308   //look through ECAL recPoints left without PRE/HCAL
309   if(ecaExist){ //if there is ecal rec point
310     Int_t iECARP ;
311     for(iECARP = 0; iECARP < nECA  ; iECARP++ ){
312       if(ecaExist[iECARP] > 0 ){
313         new ((*trackSegments)[fNTrackSegments])  
314           AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(iECARP)), null, null) ;
315         (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
316         fNTrackSegments++;    
317         if( gDebug == 2 ) 
318           Info("MakePairs", "ECAL section alone") ; 
319      } 
320     }
321   }
322   delete [] ecaExist ; 
323   delete [] preExist ; 
324   delete [] hcaExist ; 
325 }
326
327 //____________________________________________________________________________
328 void  AliEMCALTrackSegmentMakerv1::Exec(Option_t * option)
329 {
330   // STEERing method
331
332
333   if(strstr(option,"tim"))
334     gBenchmark->Start("EMCALTSMaker");
335  
336   if(strstr(option,"print")) {
337     Print("") ; 
338     return ; 
339   }
340
341   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
342
343   Int_t nevents = gime->MaxEvent() ;   
344   Int_t ievent ;
345
346   for(ievent = 0; ievent < nevents; ievent++){
347     gime->Event(ievent,"R") ;
348     //Make some initializations 
349     fNTrackSegments = 0 ;
350     
351     gime->TrackSegments()->Clear() ; 
352     
353     MakeLinks() ;
354     MakePairs() ;
355
356     WriteTrackSegments(ievent) ;
357     
358     if(strstr(option,"deb"))
359       PrintTrackSegments(option) ;
360     
361     //increment the total number of track segments per run 
362     fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ; 
363     
364   }
365   
366   if(strstr(option,"tim")){
367     gBenchmark->Stop("EMCALTSMaker");
368     Info("Exec", "took %f seconds for making TS %f seconds per event", 
369          gBenchmark->GetCpuTime("EMCALTSMaker"), gBenchmark->GetCpuTime("EMCALTSMaker")/nevents) ;
370   }
371    Unload();
372 }
373
374 //____________________________________________________________________________
375 void AliEMCALTrackSegmentMakerv1::Unload() 
376 {
377   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;  
378   gime->EmcalLoader()->UnloadRecPoints() ;
379   gime->EmcalLoader()->UnloadTracks() ;
380 }
381
382 //____________________________________________________________________________
383 void AliEMCALTrackSegmentMakerv1::Print(Option_t * option)const
384 {
385   //  Print TrackSegmentMaker parameters
386
387   Info("Print", "TrackSegmentMakerv1 parameters:") ; 
388   if( strcmp(GetName(), "") != 0 ) { 
389     printf("Making Track segments with parameters:\n") ; 
390     printf("    Allowed spred on the scalar product of two recpoints with same direction: %f\n", fClose) ;
391     printf("============================================\n") ;
392   }
393   else
394     printf("AliEMCALTrackSegmentMakerv1 not initialized ") ;
395 }
396
397 //____________________________________________________________________________
398 void AliEMCALTrackSegmentMakerv1::WriteTrackSegments(Int_t event)
399 {
400   // Writes found TrackSegments to TreeR. Creates branches 
401   // "EMCALTS" and "AliEMCALTrackSegmentMaker" with the same title.
402   // In the former branch found TrackSegments are stored, while 
403   // in the latter all parameters, with which TS were made. 
404   // ROOT does not allow overwriting existing branches, therefore
405   // first we check, if branches with the same title already exist.
406   // If yes - exits without writing.
407   
408   AliEMCALGetter *gime = AliEMCALGetter::Instance() ; 
409
410   TClonesArray * trackSegments = gime->TrackSegments() ; 
411   trackSegments->Expand(trackSegments->GetEntriesFast()) ;
412
413   TTree * treeT = gime->TreeT();
414   
415   //First TS
416   Int_t bufferSize = 32000 ;    
417   TBranch * tsBranch = treeT->Branch("EMCALTS",&trackSegments,bufferSize);
418   tsBranch->Fill() ; 
419
420   gime->WriteTracks("OVERWRITE");
421   gime->WriteTrackSegmentMaker("OVERWRITE");
422 }
423
424
425 //____________________________________________________________________________
426 void AliEMCALTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
427 {
428   // option deb - prints # of found TrackSegments
429   // option deb all - prints as well indexed of found RecParticles assigned to the TS
430   
431   TClonesArray * trackSegments = AliEMCALGetter::Instance()->TrackSegments() ; 
432
433
434   Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ; 
435   printf("nevent: %d\n", gAlice->GetEvNumber()) ; 
436   printf("        Found %d TrackSegments\n", trackSegments->GetEntriesFast() ); 
437
438   if(strstr(option,"all")) {  // printing found TS
439     printf("TrackSegment#  ECAL RP#  PRE RP#   HCAL RP#  \n") ; 
440     Int_t index;
441     for (index = 0 ; index < fNTrackSegments ; index++) {
442       AliEMCALTrackSegment * ts = (AliEMCALTrackSegment * )trackSegments->At(index) ; 
443       printf("   %d           %d        %d         %d \n", 
444              ts->GetIndexInList(), ts->GetECAIndex(), ts->GetPREIndex(), ts->GetHCAIndex() ); 
445     }   
446   }
447 }