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