new classes for track segments
[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 #include "AliRun.h"
53
54 ClassImp( AliEMCALTrackSegmentMakerv1) 
55
56
57 //____________________________________________________________________________
58   AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1() : AliEMCALTrackSegmentMaker()
59 {
60   // default ctor (to be used mainly by Streamer)
61
62   InitParameters() ; 
63
64   fTrackSegmentsInRun       = 0 ; 
65
66   fDefaultInit = kTRUE ; 
67 }
68
69 //____________________________________________________________________________
70  AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1(const char * headerFile, const char * name, const Bool_t toSplit) : AliEMCALTrackSegmentMaker(headerFile, name, toSplit)
71 {
72   // ctor
73
74   InitParameters() ; 
75
76   fTrackSegmentsInRun       = 0 ; 
77
78   Init() ;
79
80   fDefaultInit = kFALSE ; 
81
82 }
83
84 //____________________________________________________________________________
85  AliEMCALTrackSegmentMakerv1::~AliEMCALTrackSegmentMakerv1()
86
87   // dtor
88   // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
89   
90   if (!fDefaultInit) {
91     delete fPRELinkArray  ;
92     delete fHCLinkArray  ;
93     fSplitFile = 0 ; 
94   }
95 }
96
97
98 //____________________________________________________________________________
99 const TString AliEMCALTrackSegmentMakerv1::BranchName() const 
100 {  
101   TString branchName(GetName() ) ;
102   branchName.Remove(branchName.Index(Version())-1) ;
103   return branchName ;
104 }
105
106 //____________________________________________________________________________
107 Float_t  AliEMCALTrackSegmentMakerv1::HowClose(AliEMCALTowerRecPoint * ec, AliEMCALTowerRecPoint * rp, Bool_t &toofar)const
108 {
109   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
110   // Clusters are sorted in "rows" and "columns" of width 1 cm
111
112   Float_t r = -1. ;
113   Float_t delta = 10. ;  // large enough to be ineffective ??! 
114
115  
116   TVector3 vecEC = ec->XYZInAlice() ;
117   TVector3 vecRP = rp->XYZInAlice() ;
118   
119   Float_t pro = TMath::Abs(1 - (vecEC * vecRP / ( vecEC.Mag() * vecRP.Mag() ))) ; 
120
121   if(pro <= delta) { 
122     r = pro ;
123     toofar = kFALSE ;
124   } 
125   else 
126     toofar = kTRUE ;
127
128   if (gDebug == 2 ) 
129     Info("HowClose", "ec = %d, rp = %d pro = %f, toofar=%d", ec->GetIndexInList(), rp->GetIndexInList(), pro, toofar ) ; 
130  
131   return r ;
132 }
133
134 //____________________________________________________________________________
135 void  AliEMCALTrackSegmentMakerv1::Init()
136 {
137   // Make all memory allocations that are not possible in default constructor
138   
139   if ( strcmp(GetTitle(), "") == 0 )
140     SetTitle("galice.root") ;
141     
142   TString branchname = GetName() ;
143   branchname.Remove(branchname.Index(Version())-1) ;
144   AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(),branchname.Data(), fToSplit ) ; 
145   if ( gime == 0 ) {
146     Error("Init", "Could not obtain the Getter object !") ; 
147     return ;
148   } 
149   
150   fSplitFile = 0 ;
151   if(fToSplit){
152     //First - extract full path if necessary
153     TString fileName(GetTitle()) ;
154     Ssiz_t islash = fileName.Last('/') ;
155     if(islash<fileName.Length())
156       fileName.Remove(islash+1,fileName.Length()) ;
157     else
158       fileName="" ;
159     fileName+="EMCAL.RecData." ;
160     if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
161       fileName+=branchname ;
162       fileName+="." ;
163     }
164     fileName+="root" ;
165     fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));   
166     if(!fSplitFile)
167       fSplitFile =  TFile::Open(fileName.Data(),"update") ;
168   }
169   
170   fPRELinkArray = new TClonesArray("AliEMCALLink", 1000); 
171   fHCLinkArray  = new TClonesArray("AliEMCALLink", 1000); 
172
173
174   gime->PostTrackSegmentMaker(this) ;
175   gime->PostTrackSegments(BranchName()) ; 
176
177 }
178
179 //____________________________________________________________________________
180 void  AliEMCALTrackSegmentMakerv1::InitParameters()
181 {
182   fClose     = 10e-3 ;   
183   fPRELinkArray = 0 ;
184   fHCLinkArray  = 0 ;
185   TString tsmName( GetName()) ; 
186   if (tsmName.IsNull() ) 
187     tsmName = "Default" ; 
188   tsmName.Append(":") ; 
189   tsmName.Append(Version()) ; 
190   SetName(tsmName) ;
191 }
192
193
194 //____________________________________________________________________________
195 void  AliEMCALTrackSegmentMakerv1::MakeLinks()const
196
197   // Finds distances (links) between all PRE, EC and HC clusters, 
198   // which are not further apart from each other than fDangle 
199   // and sort them in accordance with this distance
200   
201   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
202   TObjArray * aECRecPoints  = gime->ECALRecPoints() ; 
203   TObjArray * aPRERecPoints = gime->PRERecPoints() ; 
204   TObjArray * aHCRecPoints  = gime->HCALRecPoints() ; 
205
206   fPRELinkArray->Clear() ;    
207   fHCLinkArray->Clear() ;    
208
209   AliEMCALTowerRecPoint * pre ;
210   AliEMCALTowerRecPoint * ec ;
211   AliEMCALTowerRecPoint * hc ;
212
213   Int_t iPRELink  = 0 ;
214   Int_t iHCLink   = 0 ;
215     
216   Int_t iECRP;
217   for(iECRP = 0; iECRP < aECRecPoints->GetEntriesFast(); iECRP++ ) {
218     ec = dynamic_cast<AliEMCALTowerRecPoint *>(aECRecPoints->At(iECRP)) ;
219     Bool_t toofar = kTRUE ;        
220     Int_t iPRERP = 0 ;    
221     for(iPRERP = 0; iPRERP < aPRERecPoints->GetEntriesFast(); iPRERP++ ) { 
222       pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(iPRERP)) ;
223       Float_t prod = HowClose(ec, pre, toofar) ;    
224       if(toofar)
225         break ;  
226       if(prod < fClose) { 
227         new ((*fPRELinkArray)[iPRELink++])  AliEMCALLink(prod, iECRP, iPRERP, 0) ;
228       }      
229     }
230     toofar = kTRUE ; 
231     Int_t iHCRP = 0 ;    
232     for(iHCRP = 0; iHCRP < aHCRecPoints->GetEntriesFast(); iHCRP++ ) { 
233       hc = dynamic_cast<AliEMCALTowerRecPoint *>(aHCRecPoints->At(iHCRP)) ;
234       Float_t prod = HowClose(ec, hc, toofar) ;    
235       if(toofar)
236         break ;  
237       if(prod < fClose) { 
238         new ((*fHCLinkArray)[iHCLink++])  AliEMCALLink(prod, iECRP, iHCRP, 1) ;
239       }      
240     }
241   } 
242   
243   fPRELinkArray->Sort() ;  //first links with largest scalar product
244   fHCLinkArray->Sort() ;  //first links with largest scalar product
245 }
246
247 //____________________________________________________________________________
248 void  AliEMCALTrackSegmentMakerv1::MakePairs()
249
250   // Using the previously made list of "links", we found the best link - i.e. 
251   // link with the largest scalar product (closest to one) to still 
252   // unassigned RecParticles. We assign these RecPoints to TrackSegment and 
253   // remove them from the list of "unassigned". 
254   
255   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
256   TObjArray * aECALRecPoints = gime->ECALRecPoints() ; 
257   TObjArray * aPRERecPoints = gime->PRERecPoints() ; 
258   TObjArray * aHCALRecPoints = gime->HCALRecPoints() ; 
259   TClonesArray * trackSegments = gime->TrackSegments(BranchName()) ;   
260     
261   //Make arrays to mark clusters already chosen
262   Int_t * ecalExist = 0;
263   Int_t nECAL = aECALRecPoints->GetEntriesFast() ;  
264   if (nECAL) 
265     ecalExist = new Int_t[nECAL] ;
266   
267   Int_t index;
268   for(index = 0; index < nECAL; index ++)
269     ecalExist[index] = 1 ;
270   
271   Bool_t * preExist = 0;
272   Int_t nPRE = aPRERecPoints->GetEntriesFast() ;  
273   if(nPRE)
274     preExist = new Bool_t[nPRE] ;
275   for(index = 0; index < nPRE; index ++)
276     preExist[index] = kTRUE ;
277   
278   Bool_t * hcalExist = 0;
279   Int_t nHCAL = aHCALRecPoints->GetEntriesFast() ;  
280   if(nHCAL)
281     hcalExist = new Bool_t[nHCAL] ;
282   for(index = 0; index < nHCAL; index ++)
283     hcalExist[index] = kTRUE ;
284   
285   AliEMCALTowerRecPoint * null = 0 ; 
286  // Finds the smallest links and makes pairs of PRE and ECAL clusters with largest scalar product 
287  
288   TIter nextPRE(fPRELinkArray) ;
289   AliEMCALLink * linkPRE ;
290   
291   while ( (linkPRE =  static_cast<AliEMCALLink *>(nextPRE()) ) ){  
292
293     if(ecalExist[linkPRE->GetECAL()] != -1){ //without PRE yet 
294
295       if(preExist[linkPRE->GetOther()]){ // PRE still exist
296         
297         new ((* trackSegments)[fNTrackSegments]) 
298           AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECALRecPoints->At(linkPRE->GetECAL())) , 
299                                dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(linkPRE->GetOther())), null) ;
300         (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
301         fNTrackSegments++ ;
302         if (gDebug == 2 ) 
303           Info("MakePairs", "ECAL section with PRE section") ;  
304         ecalExist[linkPRE->GetECAL()] = -1 ; //Mark ecal  that pre was found 
305         //mark PRE recpoint as already used 
306         preExist[linkPRE->GetOther()] = kFALSE ;
307       } //if PRE still exist
308     } 
309   }
310
311   // Finds the smallest links and makes pairs of HCAL and ECAL clusters with largest scalar product 
312  
313   TIter nextHCAL(fHCLinkArray) ;
314   AliEMCALLink * linkHCAL ;
315   
316   while ( (linkHCAL =  static_cast<AliEMCALLink *>(nextHCAL()) ) ){  
317
318     if(ecalExist[linkHCAL->GetECAL()] != -2){ //without HCAL yet 
319
320       if(hcalExist[linkHCAL->GetOther()]){ // HCAL still exist
321         // search among the already existing track segments
322         Int_t ii ; 
323         Bool_t found = kFALSE ; 
324         AliEMCALTrackSegment * ts = 0 ; 
325         for ( ii = 0 ; ii < fNTrackSegments ; ii++ ) {
326           ts = dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(ii)) ;
327           if ( ts->GetECIndex() == linkHCAL->GetECAL() ) {
328             found = kTRUE ; 
329             break ; 
330           }
331         }
332         if (found){
333           ts->SetHCRecPoint( dynamic_cast<AliEMCALTowerRecPoint *>(aHCALRecPoints->At(linkHCAL->GetOther())) ) ;
334           if (gDebug == 2 ) 
335             Info("MakePairs", "ECAL section with PRE and HCAL sections") ;
336         }       
337         if (!found) {
338           new ((* trackSegments)[fNTrackSegments]) 
339             AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECALRecPoints->At(linkHCAL->GetECAL())), null,
340                                  dynamic_cast<AliEMCALTowerRecPoint *>(aHCALRecPoints->At(linkHCAL->GetOther()))) ; 
341         (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
342         fNTrackSegments++ ;
343         if (gDebug == 2 ) 
344           Info("MakePairs", "ECAL section with HCAL section") ;         
345         }
346         ecalExist[linkHCAL->GetECAL()] = -2 ; //Mark ecal  that hcal was found 
347         //mark HCAL recpoint as already used 
348         hcalExist[linkHCAL->GetOther()] = kFALSE ;
349       } //if HCAL still exist
350     } 
351   }
352   
353
354   //look through ECAL recPoints left without PRE/HCAL
355   if(ecalExist){ //if there is ecal rec point
356     Int_t iECALRP ;
357     for(iECALRP = 0; iECALRP < nECAL  ; iECALRP++ ){
358       if(ecalExist[iECALRP] > 0 ){
359         new ((*trackSegments)[fNTrackSegments])  
360           AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECALRecPoints->At(iECALRP)), null, null) ;
361         (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
362         fNTrackSegments++;    
363         if( gDebug == 2 ) 
364           Info("MakePairs", "ECAL section alone") ; 
365      } 
366     }
367   }
368   delete [] ecalExist ; 
369   delete [] preExist ; 
370   delete [] hcalExist ; 
371 }
372
373 //____________________________________________________________________________
374 void  AliEMCALTrackSegmentMakerv1::Exec(Option_t * option)
375 {
376   // STEERing method
377
378   if( strcmp(GetName(), "")== 0 ) 
379     Init() ;
380
381   if(strstr(option,"tim"))
382     gBenchmark->Start("EMCALTSMaker");
383  
384   if(strstr(option,"print")) {
385     Print("") ; 
386     return ; 
387   }
388
389   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
390   if(gime->BranchExists("TrackSegments") )
391     return ;
392
393   Int_t nevents = gime->MaxEvent() ;       //(Int_t) gAlice->TreeE()->GetEntries() ;
394   Int_t ievent ;
395
396   for(ievent = 0; ievent < nevents; ievent++){
397
398     gime->Event(ievent,"R") ;
399  
400     //Make some initializations 
401     fNTrackSegments = 0 ;
402     gime->TrackSegments(BranchName())->Clear() ; 
403     
404     MakeLinks() ;
405     
406     MakePairs() ;
407
408     WriteTrackSegments(ievent) ;
409     
410     if(strstr(option,"deb"))
411       PrintTrackSegments(option) ;
412     
413     //increment the total number of track segments per run 
414     fTrackSegmentsInRun += gime->TrackSegments(BranchName())->GetEntriesFast() ; 
415     
416   }
417   
418   if(strstr(option,"tim")){
419     gBenchmark->Stop("EMCALTSMaker");
420     Info("Exec", "took %f seconds for making TS %f seconds per event", 
421          gBenchmark->GetCpuTime("EMCALTSMaker"), gBenchmark->GetCpuTime("EMCALTSMaker")/nevents) ;
422   }
423   
424 }
425
426 //____________________________________________________________________________
427 void AliEMCALTrackSegmentMakerv1::Print(Option_t * option)const
428 {
429   //  Print TrackSegmentMaker parameters
430
431   Info("Print", "TrackSegmentMakerv1 parameters:") ; 
432   if( strcmp(GetName(), "") != 0 ) { 
433     printf("Making Track segments with parameters:\n") ; 
434     printf("    Allowed spred on the scalar product of two recpoints with same direction: %f\n", fClose) ;
435     printf("============================================\n") ;
436   }
437   else
438     printf("AliEMCALTrackSegmentMakerv1 not initialized ") ;
439 }
440
441 //____________________________________________________________________________
442 void AliEMCALTrackSegmentMakerv1::WriteTrackSegments(Int_t event)
443 {
444   // Writes found TrackSegments to TreeR. Creates branches 
445   // "EMCALTS" and "AliEMCALTrackSegmentMaker" with the same title.
446   // In the former branch found TrackSegments are stored, while 
447   // in the latter all parameters, with which TS were made. 
448   // ROOT does not allow overwriting existing branches, therefore
449   // first we check, if branches with the same title already exist.
450   // If yes - exits without writing.
451   
452   AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; 
453
454   TClonesArray * trackSegments = gime->TrackSegments() ; 
455   trackSegments->Expand(trackSegments->GetEntriesFast()) ;
456   TTree * treeR ;
457   
458   if(fToSplit){
459     if(!fSplitFile)
460       return ;
461     fSplitFile->cd() ;
462     char name[10] ;
463     sprintf(name,"%s%d", "TreeR",event) ;
464     treeR = dynamic_cast<TTree*>(fSplitFile->Get(name)); 
465   }
466   else{
467     treeR = gAlice->TreeR();
468   }
469   
470   if(!treeR){
471     gAlice->MakeTree("R", fSplitFile);
472     treeR = gAlice->TreeR() ;
473   }
474   
475   //First TS
476   Int_t bufferSize = 32000 ;    
477   TBranch * tsBranch = treeR->Branch("EMCALTS",&trackSegments,bufferSize);
478   tsBranch->SetTitle(BranchName());
479
480   //Second -TSMaker
481   Int_t splitlevel = 0 ;
482   AliEMCALTrackSegmentMakerv1 * ts = this ;
483   TBranch * tsMakerBranch = treeR->Branch("AliEMCALTrackSegmentMaker","AliEMCALTrackSegmentMakerv1",
484                                           &ts,bufferSize,splitlevel);
485   tsMakerBranch->SetTitle(BranchName());
486
487   tsBranch->Fill() ;  
488   tsMakerBranch->Fill() ;
489
490   treeR->AutoSave() ; //Write(0,kOverwrite) ;  
491   if(gAlice->TreeR()!=treeR)
492     treeR->Delete();
493 }
494
495
496 //____________________________________________________________________________
497 void AliEMCALTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
498 {
499   // option deb - prints # of found TrackSegments
500   // option deb all - prints as well indexed of found RecParticles assigned to the TS
501   TString taskName(GetName()) ; 
502   taskName.Remove(taskName.Index(Version())-1) ;
503   
504   TClonesArray * trackSegments = AliEMCALGetter::GetInstance()->TrackSegments(taskName) ; 
505
506
507   Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ; 
508   printf("nevent: %d\n", gAlice->GetEvNumber()) ; 
509   printf("        Found %d TrackSegments\n", trackSegments->GetEntriesFast() ); 
510
511   if(strstr(option,"all")) {  // printing found TS
512     printf("TrackSegment#  ECAL RP#  PRE RP#   HCAL RP#  \n") ; 
513     Int_t index;
514     for (index = 0 ; index < fNTrackSegments ; index++) {
515       AliEMCALTrackSegment * ts = (AliEMCALTrackSegment * )trackSegments->At(index) ; 
516       printf("   %d           %d        %d         %d \n", 
517              ts->GetIndexInList(), ts->GetECIndex(), ts->GetPREIndex(), ts->GetHCIndex() ); 
518     }   
519   }
520 }