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() ;
 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()
 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 }