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 }