PHOS/AliPHOSTrackSegmentMakerv1.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 PHOS track segments
  18 // Track segment for PHOS is list of
  19 //        EMC RecPoint + (possibly) CPV RecPoint + (possibly) PPSD RecPoint
  20 // To find TrackSegments we do the following:
  21 //  for each EMC RecPoints we look at
  22 //   CPV/PPSD RecPoints in the radious fRcpv.
  23 //  If there is such a CPV RecPoint,
  24 //   we make "Link" it is just indexes of EMC and CPV/PPSD RecPoint and distance
  25 //   between them in the PHOS plane.
  26 //  Then we sort "Links" and starting from the
  27 //   least "Link" pointing to the unassined EMC and CPV RecPoints assing them to
  28 //   new TrackSegment.
  29 // If there is no CPV/PPSD RecPoint we make TrackSegment
  30 // consisting from EMC alone. There is no TrackSegments without EMC RecPoint.
  31 //// In principle this class should be called from AliPHOSReconstructioner, but
  32 // one can use it as well in standalone mode.
  33 // Use  case:
  34 //  root [0] AliPHOSTrackSegmentMakerv1 * t = new AliPHOSTrackSegmentMaker("galice.root", "tracksegmentsname", "recpointsname")
  35 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  36 //               // reads gAlice from header file "galice.root", uses recpoints stored in the branch names "recpointsname" (default = "Default")
  37 //               // and saves recpoints in branch named "tracksegmentsname" (default = "recpointsname")
  38 //  root [1] t->ExecuteTask()
  39 //  root [2] t->SetMaxEmcPpsdDistance(5)
  40 //  root [3] t->SetTrackSegmentsBranch("max distance 5 cm")
  41 //  root [4] t->ExecuteTask("deb all time")
  42 //
  43 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH)
  44 //
  45
  46 // --- ROOT system ---
  47 #include "TTree.h"
  48 #include "TBenchmark.h"
  49
  50 // --- Standard library ---
  51
  52 // --- AliRoot header files ---
  53 #include "AliPHOSGeometry.h"
  54 #include "AliPHOSTrackSegmentMakerv1.h"
  55 #include "AliPHOSTrackSegment.h"
  56 #include "AliPHOSLink.h"
  57 #include "AliPHOSGetter.h"
  58
  59 ClassImp( AliPHOSTrackSegmentMakerv1)
  60
  61
  62 //____________________________________________________________________________
  63   AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() : AliPHOSTrackSegmentMaker()
  64 {
  65   // default ctor (to be used mainly by Streamer)
  66
  67   InitParameters() ;
  68   fDefaultInit = kTRUE ;
  69 }
  70
  71 //____________________________________________________________________________
  72  AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const TString alirunFileName, const TString eventFolderName)
  73    :AliPHOSTrackSegmentMaker(alirunFileName, eventFolderName)
  74 {
  75   // ctor
  76
  77   InitParameters() ;
  78   Init() ;
  79   fDefaultInit = kFALSE ;
  80 }
  81
  82 //____________________________________________________________________________
  83  AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
  84 {
  85   // dtor
  86   // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
  87   if (!fDefaultInit)
  88     delete fLinkUpArray ;
  89 }
  90
  91
  92 //____________________________________________________________________________
  93 const TString AliPHOSTrackSegmentMakerv1::BranchName() const
  94 {
  95
  96   return GetName() ;
  97 }
  98
  99 //____________________________________________________________________________
 100 void  AliPHOSTrackSegmentMakerv1::FillOneModule()
 101 {
 102   // Finds first and last indexes between which
 103   // clusters from one PHOS module are
 104
 105   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 106
 107   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 108   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 109
 110   //First EMC clusters
 111   Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
 112   for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
 113         ((dynamic_cast<AliPHOSRecPoint *>(emcRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule );
 114       fEmcLast ++)  ;
 115
 116   //Now CPV clusters
 117   Int_t totalCpv = cpvRecPoints->GetEntriesFast() ;
 118
 119     for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) &&
 120          ((dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule );
 121        fCpvLast ++) ;
 122
 123 }
 124
 125 //____________________________________________________________________________
 126 Float_t  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,AliPHOSRecPoint * cpvClu, Bool_t &toofar)const
 127 {
 128   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
 129   // Clusters are sorted in "rows" and "columns" of width 1 cm
 130
 131   Float_t delta = 1 ;  // Width of the rows in sorting of RecPoints (in cm)
 132                        // if you change this value, change it as well in xxxRecPoint::Compare()
 133   Float_t r = fRcpv ;
 134
 135   TVector3 vecEmc ;
 136   TVector3 vecCpv ;
 137
 138   emcClu->GetLocalPosition(vecEmc) ;
 139   cpvClu->GetLocalPosition(vecCpv)  ;
 140
 141   if(emcClu->GetPHOSMod() == cpvClu->GetPHOSMod()){
 142     if(vecCpv.X() <= vecEmc.X() + fRcpv + 2*delta ){
 143
 144       vecCpv = vecCpv  - vecEmc ;
 145       r = vecCpv.Mag() ;
 146       toofar = kFALSE ;
 147
 148     } // if  xPpsd >= xEmc + ...
 149     else
 150       toofar = kTRUE ;
 151   }
 152   else
 153     toofar = kTRUE ;
 154
 155   //toofar = kFALSE ;
 156
 157
 158   return r ;
 159 }
 160
 161 //____________________________________________________________________________
 162 void  AliPHOSTrackSegmentMakerv1::Init()
 163 {
 164   // Make all memory allocations that are not possible in default constructor
 165
 166   AliPHOSGetter* gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
 167
 168   fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000);
 169   if ( !gime->TrackSegmentMaker() ) {
 170     gime->PostTrackSegmentMaker(this);
 171   }
 172 }
 173
 174 //____________________________________________________________________________
 175 void  AliPHOSTrackSegmentMakerv1::InitParameters()
 176 {
 177   //Initializes parameters
 178   fRcpv      = 10. ;
 179   fEmcFirst  = 0 ;
 180   fEmcLast   = 0 ;
 181   fCpvFirst  = 0 ;
 182   fCpvLast   = 0 ;
 183   fLinkUpArray = 0 ;
 184   fTrackSegmentsInRun       = 0 ;
 185 }
 186
 187
 188 //____________________________________________________________________________
 189 void  AliPHOSTrackSegmentMakerv1::MakeLinks()const
 190 {
 191   // Finds distances (links) between all EMC and PPSD clusters,
 192   // which are not further apart from each other than fRcpv
 193   // and sort them in accordance with this distance
 194
 195   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 196   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 197   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 198
 199   fLinkUpArray->Clear() ;
 200
 201   AliPHOSRecPoint * cpv ;
 202   AliPHOSEmcRecPoint * emcclu ;
 203
 204   Int_t iLinkUp  = 0 ;
 205
 206   Int_t iEmcRP;
 207   for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
 208     emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP)) ;
 209
 210     Bool_t toofar ;
 211     Int_t iCpv = 0 ;
 212     for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) {
 213
 214       cpv = dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(iCpv)) ;
 215       Float_t r = GetDistanceInPHOSPlane(emcclu, cpv, toofar) ;
 216
 217       if(toofar)
 218         break ;
 219       if(r < fRcpv) {
 220         new ((*fLinkUpArray)[iLinkUp++])  AliPHOSLink(r, iEmcRP, iCpv) ;
 221       }
 222     }
 223   }
 224
 225   fLinkUpArray->Sort() ;  //first links with smallest distances
 226 }
 227
 228 //____________________________________________________________________________
 229 void  AliPHOSTrackSegmentMakerv1::MakePairs()
 230 {
 231   // Using the previously made list of "links", we found the smallest link - i.e.
 232   // link with the least distance between EMC and CPV and pointing to still
 233   // unassigned RecParticles. We assign these RecPoints to TrackSegment and
 234   // remove them from the list of "unassigned".
 235
 236   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 237
 238   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 239   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 240   TClonesArray * trackSegments = gime->TrackSegments();
 241
 242   //Make arrays to mark clusters already chosen
 243   Int_t * emcExist = 0;
 244   if(fEmcLast > fEmcFirst)
 245     emcExist = new Int_t[fEmcLast-fEmcFirst] ;
 246
 247   Int_t index;
 248   for(index = 0; index <fEmcLast-fEmcFirst; index ++)
 249     emcExist[index] = 1 ;
 250
 251   Bool_t * cpvExist = 0;
 252   if(fCpvLast > fCpvFirst)
 253     cpvExist = new Bool_t[fCpvLast-fCpvFirst] ;
 254   for(index = 0; index <fCpvLast-fCpvFirst; index ++)
 255     cpvExist[index] = kTRUE ;
 256
 257
 258   // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance
 259   TIter nextUp(fLinkUpArray) ;
 260
 261   AliPHOSLink * linkUp ;
 262
 263   AliPHOSRecPoint * nullpointer = 0 ;
 264
 265   while ( (linkUp =  static_cast<AliPHOSLink *>(nextUp()) ) ){
 266
 267     if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){ //without ppsd Up yet
 268
 269       if(cpvExist[linkUp->GetPpsd()-fCpvFirst]){ //CPV still exist
 270
 271        new ((* trackSegments)[fNTrackSegments])
 272          AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(linkUp->GetEmc())) ,
 273                            dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(linkUp->GetPpsd()))) ;
 274        (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
 275        fNTrackSegments++ ;
 276
 277        emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc  that Cpv was found
 278        //mark CPV recpoint as already used
 279        cpvExist[linkUp->GetPpsd()-fCpvFirst] = kFALSE ;
 280       } //if ppsdUp still exist
 281     }
 282   }
 283
 284   //look through emc recPoints left without CPV/PPSD
 285   if(emcExist){ //if there is emc rec point
 286     Int_t iEmcRP ;
 287     for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst  ; iEmcRP++ ){
 288       if(emcExist[iEmcRP] > 0 ){
 289        new ((*trackSegments)[fNTrackSegments])
 290          AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP+fEmcFirst)),
 291                            nullpointer) ;
 292        (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
 293        fNTrackSegments++;
 294       }
 295     }
 296   }
 297   delete [] emcExist ;
 298   delete [] cpvExist ;
 299 }
 300
 301 //____________________________________________________________________________
 302 void  AliPHOSTrackSegmentMakerv1::Exec(Option_t * option)
 303 {
 304   // STEERing method
 305
 306   if(strstr(option,"tim"))
 307     gBenchmark->Start("PHOSTSMaker");
 308
 309   if(strstr(option,"print")) {
 310     Print() ;
 311     return ;
 312   }
 313
 314   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 315
 316   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 317
 318   Int_t nevents = gime->MaxEvent() ;
 319   Int_t ievent ;
 320
 321   for(ievent = 0; ievent < nevents; ievent++) {
 322     gime->Event(ievent,"R") ;
 323     //Make some initializations
 324     fNTrackSegments = 0 ;
 325     fEmcFirst = 0 ;
 326     fEmcLast  = 0 ;
 327     fCpvFirst = 0 ;
 328     fCpvLast  = 0 ;
 329
 330     gime->TrackSegments()->Clear();
 331
 332     //    if(!ReadRecPoints(ievent))   continue; //reads RecPoints for event ievent
 333
 334     for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ) {
 335       FillOneModule() ;
 336       MakeLinks() ;
 337       MakePairs() ;
 338     }
 339
 340     WriteTrackSegments() ;
 341
 342     if(strstr(option,"deb"))
 343       PrintTrackSegments(option);
 344
 345     //increment the total number of track segments per run
 346     fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ;
 347
 348   }
 349
 350   if(strstr(option,"tim")){
 351     gBenchmark->Stop("PHOSTSMaker");
 352     Info("Exec", "took %f seconds for making TS %f seconds per event",
 353           gBenchmark->GetCpuTime("PHOSTSMaker"),
 354           gBenchmark->GetCpuTime("PHOSTSMaker")/nevents) ;
 355    }
 356   Unload();
 357 }
 358
 359 //____________________________________________________________________________
 360 void AliPHOSTrackSegmentMakerv1::Unload()
 361 {
 362   // Unloads the task from the folder
 363   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 364   gime->PhosLoader()->UnloadRecPoints() ;
 365   gime->PhosLoader()->UnloadTracks() ;
 366 }
 367
 368 //____________________________________________________________________________
 369 void AliPHOSTrackSegmentMakerv1::Print()const
 370 {
 371   //  Print TrackSegmentMaker parameters
 372
 373   TString message("") ;
 374   if( strcmp(GetName(), "") != 0 ) {
 375     message = "\n======== AliPHOSTrackSegmentMakerv1 ========\n" ;
 376     message += "Making Track segments\n" ;
 377     message += "with parameters:\n" ;
 378     message += "     Maximal EMC - CPV (PPSD) distance (cm) %f\n" ;
 379     message += "============================================\n" ;
 380     Info("Print", message.Data(),fRcpv) ;
 381   }
 382   else
 383     Info("Print", "AliPHOSTrackSegmentMakerv1 not initialized ") ;
 384 }
 385
 386 //____________________________________________________________________________
 387 void AliPHOSTrackSegmentMakerv1::WriteTrackSegments()
 388 {
 389   // Writes found TrackSegments to TreeR. Creates branches
 390   // "PHOSTS" and "AliPHOSTrackSegmentMaker" 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   AliPHOSGetter *gime = AliPHOSGetter::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("PHOSTS",&trackSegments,bufferSize);
 407   tsBranch->Fill() ;
 408
 409   gime->WriteTracks("OVERWRITE");
 410   gime->WriteTrackSegmentMaker("OVERWRITE");
 411 }
 412
 413
 414 //____________________________________________________________________________
 415 void AliPHOSTrackSegmentMakerv1::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 = AliPHOSGetter::Instance()->TrackSegments() ;
 421
 422   Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
 423   printf("nevent: %d\n", gAlice->GetEvNumber()) ;
 424   printf("        Found %d TrackSegments\n", trackSegments->GetEntriesFast() );
 425
 426   if(strstr(option,"all")) {  // printing found TS
 427     printf("TrackSegment #  EMC RP#  CPV RP#\n") ;
 428     Int_t index;
 429     for (index = 0 ; index <trackSegments->GetEntriesFast() ; index++) {
 430       AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )trackSegments->At(index) ;
 431       printf("   %d           %d        %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetCpvIndex() ) ;
 432     }
 433   }
 434 }