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 }