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(AliEMCALTowerRecPoint * ec, AliEMCALTowerRecPoint * 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 = ec->XYZInAlice() ;
  97   TVector3 vecRP = rp->XYZInAlice() ;
  98
  99   Float_t pro = TMath::Abs(1 - (vecEC * vecRP / ( vecEC.Mag() * vecRP.Mag() ))) ;
 100
 101   if(pro <= delta) {
 102     r = pro ;
 103     toofar = kFALSE ;
 104   }
 105   else
 106     toofar = kTRUE ;
 107
 108   if (gDebug == 2 )
 109     printf("HowClose: ec = %d, rp = %d pro = %f, toofar=%d", ec->GetIndexInList(), rp->GetIndexInList(), pro, toofar ) ;
 110
 111   return r ;
 112 }
 113
 114 //____________________________________________________________________________
 115 void  AliEMCALTrackSegmentMakerv1::Init()
 116 {
 117   // Make all memory allocations that are not possible in default constructor
 118
 119   AliEMCALGetter* gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());
 120
 121   if ( !gime->TrackSegmentMaker() ) {
 122     gime->PostTrackSegmentMaker(this);
 123   }
 124 }
 125
 126 //____________________________________________________________________________
 127 void  AliEMCALTrackSegmentMakerv1::InitParameters()
 128 {
 129   fClose              = 10e-3 ;
 130   fTrackSegmentsInRun = 0 ;
 131 }
 132
 133
 134 //____________________________________________________________________________
 135 void  AliEMCALTrackSegmentMakerv1::MakeLinks()const
 136 {
 137   // Finds distances (links) between all PRE, EC and HC clusters,
 138   // which are not further apart from each other than fDangle
 139   // and sort them in accordance with this distance
 140
 141   /* AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 142   TObjArray * aECARecPoints  = gime->ECARecPoints() ;
 143   // TObjArray * aPRERecPoints  = gime->PRERecPoints() ;
 144   //TObjArray * aHCARecPoints  = gime->HCARecPoints() ;
 145
 146   fPRELinkArray->Clear() ;
 147   fHCALinkArray->Clear() ;
 148
 149   AliEMCALTowerRecPoint * pre ;
 150   AliEMCALTowerRecPoint * eca ;
 151   AliEMCALTowerRecPoint * hca ;
 152
 153   Int_t iPRELink  = 0 ;
 154   Int_t iHCALink  = 0 ;
 155
 156   Int_t iECARP;
 157   for(iECARP = 0; iECARP < aECARecPoints->GetEntriesFast(); iECARP++ ) {
 158     eca = dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(iECARP)) ;
 159     Bool_t toofar = kTRUE ;
 160     Int_t iPRERP = 0 ;
 161     for(iPRERP = 0; iPRERP < aPRERecPoints->GetEntriesFast(); iPRERP++ ) {
 162       pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(iPRERP)) ;
 163       Float_t prod = HowClose(eca, pre, toofar) ;
 164       if(toofar)
 165         break ;
 166       if(prod < fClose) {
 167         new ((*fPRELinkArray)[iPRELink++])  AliEMCALLink(prod, iECARP, iPRERP, 0) ;
 168       }
 169     }
 170     toofar = kTRUE ;
 171     Int_t iHCARP = 0 ;
 172     for(iHCARP = 0; iHCARP < aHCARecPoints->GetEntriesFast(); iHCARP++ ) {
 173       hca = dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(iHCARP)) ;
 174       Float_t prod = HowClose(eca, hca, toofar) ;
 175       if(toofar)
 176         break ;
 177       if(prod < fClose) {
 178         new ((*fHCALinkArray)[iHCALink++])  AliEMCALLink(prod, iECARP, iHCARP, 1) ;
 179       }
 180     }
 181   }
 182
 183   fPRELinkArray->Sort() ;  //first links with largest scalar product
 184   fHCALinkArray->Sort() ;  //first links with largest scalar product
 185   */
 186 }
 187
 188 //____________________________________________________________________________
 189 void  AliEMCALTrackSegmentMakerv1::MakePairs()
 190 {
 191   // Using the previously made list of "links", we found the best link - i.e.
 192   // link with the largest scalar product (closest to one) to still
 193   // unassigned RecParticles. We assign these RecPoints to TrackSegment and
 194   // remove them from the list of "unassigned".
 195
 196   /*AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 197    TObjArray * aECARecPoints = gime->ECARecPoints() ;
 198   TObjArray * aPRERecPoints = gime->PRERecPoints() ;
 199   TObjArray * aHCARecPoints = gime->HCARecPoints() ;
 200   TClonesArray * trackSegments = gime->TrackSegments() ;
 201
 202   //Make arrays to mark clusters already chosen
 203   Int_t * ecaExist = 0;
 204   Int_t nECA = aECARecPoints->GetEntriesFast() ;
 205   if (nECA)
 206     ecaExist = new Int_t[nECA] ;
 207
 208   Int_t index;
 209   for(index = 0; index < nECA; index ++)
 210     ecaExist[index] = 1 ;
 211
 212   Bool_t * preExist = 0;
 213   Int_t nPRE = aPRERecPoints->GetEntriesFast() ;
 214   if(nPRE)
 215     preExist = new Bool_t[nPRE] ;
 216   for(index = 0; index < nPRE; index ++)
 217     preExist[index] = kTRUE ;
 218
 219   Bool_t * hcaExist = 0;
 220   Int_t nHCA = aHCARecPoints->GetEntriesFast() ;
 221   if(nHCA)
 222     hcaExist = new Bool_t[nHCA] ;
 223   for(index = 0; index < nHCA; index ++)
 224     hcaExist[index] = kTRUE ;
 225
 226   AliEMCALTowerRecPoint * null = 0 ;
 227  // Finds the smallest links and makes pairs of PRE and ECAL clusters with largest scalar product
 228
 229   TIter nextPRE(fPRELinkArray) ;
 230   AliEMCALLink * linkPRE ;
 231
 232   while ( (linkPRE =  static_cast<AliEMCALLink *>(nextPRE()) ) ){
 233
 234     if(ecaExist[linkPRE->GetECA()] != -1){ //without PRE yet
 235
 236       if(preExist[linkPRE->GetOther()]){ // PRE still exist
 237
 238         new ((* trackSegments)[fNTrackSegments])
 239           AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(linkPRE->GetECA())) ,
 240                                dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(linkPRE->GetOther())), null) ;
 241         (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
 242         fNTrackSegments++ ;
 243         if (gDebug == 2 )
 244           printf("MakePairs: ECAL section with PRE section") ;
 245         ecaExist[linkPRE->GetECA()] = -1 ; //Mark ecal  that pre was found
 246         //mark PRE recpoint as already used
 247         preExist[linkPRE->GetOther()] = kFALSE ;
 248       } //if PRE still exist
 249     }
 250   }
 251
 252   // Finds the smallest links and makes pairs of HCAL and ECAL clusters with largest scalar product
 253
 254   TIter nextHCA(fHCALinkArray) ;
 255   AliEMCALLink * linkHCA ;
 256
 257   while ( (linkHCA =  static_cast<AliEMCALLink *>(nextHCA()) ) ){
 258
 259     if(ecaExist[linkHCA->GetECA()] != -2){ //without HCAL yet
 260
 261       if(hcaExist[linkHCA->GetOther()]){ // HCAL still exist
 262         // search among the already existing track segments
 263         Int_t ii ;
 264         Bool_t found = kFALSE ;
 265         AliEMCALTrackSegment * ts = 0 ;
 266         for ( ii = 0 ; ii < fNTrackSegments ; ii++ ) {
 267           ts = dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(ii)) ;
 268           if ( ts->GetECAIndex() == linkHCA->GetECA() ) {
 269             found = kTRUE ;
 270             break ;
 271           }
 272         }
 273         if (found){
 274           ts->SetHCARecPoint( dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(linkHCA->GetOther())) ) ;
 275           if (gDebug == 2 )
 276             printf("MakePairs: ECAL section with PRE and HCAL sections") ;
 277         }
 278         if (!found) {
 279           new ((* trackSegments)[fNTrackSegments])
 280             AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(linkHCA->GetECA())), null,
 281                                  dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(linkHCA->GetOther()))) ;
 282         (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
 283         fNTrackSegments++ ;
 284         if (gDebug == 2 )
 285           printf("MakePairs: ECAL section with HCAL section") ;
 286         }
 287         ecaExist[linkHCA->GetECA()] = -2 ; //Mark ecal  that hcal was found
 288         //mark HCAL recpoint as already used
 289         hcaExist[linkHCA->GetOther()] = kFALSE ;
 290       } //if HCAL still exist
 291     }
 292   }
 293
 294
 295   //look through ECAL recPoints left without PRE/HCAL
 296   if(ecaExist){ //if there is ecal rec point
 297     Int_t iECARP ;
 298     for(iECARP = 0; iECARP < nECA  ; iECARP++ ){
 299       if(ecaExist[iECARP] > 0 ){
 300         new ((*trackSegments)[fNTrackSegments])
 301           AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(iECARP)), null, null) ;
 302         (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
 303         fNTrackSegments++;
 304         if( gDebug == 2 )
 305           printf("MakePairs: ECAL section alone") ;
 306      }
 307     }
 308   }
 309   delete [] ecaExist ;
 310   delete [] preExist ;
 311   delete [] hcaExist ;
 312   */
 313 }
 314
 315 //____________________________________________________________________________
 316 void  AliEMCALTrackSegmentMakerv1::Exec(Option_t * option)
 317 {
 318   // STEERing method
 319
 320
 321   if(strstr(option,"tim"))
 322     gBenchmark->Start("EMCALTSMaker");
 323
 324   if(strstr(option,"print")) {
 325     Print("") ;
 326     return ;
 327   }
 328
 329   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 330
 331   Int_t nevents = gime->MaxEvent() ;
 332   Int_t ievent ;
 333
 334   for(ievent = 0; ievent < nevents; ievent++){
 335     gime->Event(ievent,"R") ;
 336     //Make some initializations
 337     fNTrackSegments = 0 ;
 338
 339     gime->TrackSegments()->Clear() ;
 340
 341     MakeLinks() ;
 342     MakePairs() ;
 343
 344     WriteTrackSegments() ;
 345
 346     if(strstr(option,"deb"))
 347       PrintTrackSegments(option) ;
 348
 349     //increment the total number of track segments per run
 350     fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ;
 351
 352   }
 353
 354   if(strstr(option,"tim")){
 355     gBenchmark->Stop("EMCALTSMaker");
 356     printf("Exec: took %f seconds for making TS %f seconds per event",
 357          gBenchmark->GetCpuTime("EMCALTSMaker"), gBenchmark->GetCpuTime("EMCALTSMaker")/nevents) ;
 358   }
 359    Unload();
 360 }
 361
 362 //____________________________________________________________________________
 363 void AliEMCALTrackSegmentMakerv1::Unload()
 364 {
 365   // Unloads the RecPoints and Tracks
 366   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 367   gime->EmcalLoader()->UnloadRecPoints() ;
 368   gime->EmcalLoader()->UnloadTracks() ;
 369 }
 370
 371 //____________________________________________________________________________
 372 void AliEMCALTrackSegmentMakerv1::Print(Option_t * /*option*/)const
 373 {
 374   //  Print TrackSegmentMaker parameters
 375
 376   printf("Print: TrackSegmentMakerv1 parameters:") ;
 377   if( strcmp(GetName(), "") != 0 ) {
 378     printf("Making Track segments with parameters:\n") ;
 379     printf("    Allowed spred on the scalar product of two recpoints with same direction: %f\n", fClose) ;
 380     printf("============================================\n") ;
 381   }
 382   else
 383     printf("AliEMCALTrackSegmentMakerv1 not initialized ") ;
 384 }
 385
 386 //____________________________________________________________________________
 387 void AliEMCALTrackSegmentMakerv1::WriteTrackSegments()
 388 {
 389   // Writes found TrackSegments to TreeR. Creates branches
 390   // "EMCALTS" and "AliEMCALTrackSegmentMaker" 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   AliEMCALGetter *gime = AliEMCALGetter::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("EMCALTS",&trackSegments,bufferSize);
 407   tsBranch->Fill() ;
 408
 409   gime->WriteTracks("OVERWRITE");
 410   gime->WriteTrackSegmentMaker("OVERWRITE");
 411 }
 412
 413
 414 //____________________________________________________________________________
 415 void AliEMCALTrackSegmentMakerv1::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 = AliEMCALGetter::Instance()->TrackSegments() ;
 421
 422
 423   printf("PrintTrackSegments: Results from TrackSegmentMaker:") ;
 424   printf("nevent: %d\n", gAlice->GetEvNumber()) ;
 425   printf("        Found %d TrackSegments\n", trackSegments->GetEntriesFast() );
 426
 427   if(strstr(option,"all")) {  // printing found TS
 428     printf("TrackSegment#  ECAL RP#  \n") ;
 429     Int_t index;
 430     for (index = 0 ; index < fNTrackSegments ; index++) {
 431       AliEMCALTrackSegment * ts = (AliEMCALTrackSegment * )trackSegments->At(index) ;
 432       printf("   %d           %d \n",
 433              ts->GetIndexInList(), ts->GetECAIndex());
 434     }
 435   }
 436 }