PHOS/AliPHOSPIDv0.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
  16 /* $Id$ */
  17
  18 /* History of cvs commits:
  19  *
  20  * $Log$
  21  */
  22
  23 //_________________________________________________________________________
  24 // Implementation version v0 of the PHOS particle identifier
  25 // Particle identification based on the
  26 //     - CPV information,
  27 //     - Preshower information (in MIXT or GPS2 geometries)
  28 //     - shower width.
  29 //
  30 // CPV or Preshower clusters should be closer in PHOS plane than fCpvEmcDistance (in cm).
  31 // This parameter can be set by method SetCpvtoEmcDistanceCut(Float_t cut)
  32 //
  33 // One can set desirable ID method by the function SetIdentificationMethod(option).
  34 // Presently the following options can be used together or separately :
  35 //     - "disp": use dispersion cut on shower width
  36 //               (width can be set by method SetDispersionCut(Float_t cut)
  37 //     - "ell" : use cut on the axis of the ellipse, drawn around shower
  38 //       (this cut can be changed by SetShowerProfileCut(char* formula),
  39 //        where formula - any function of two variables f(lambda[0],lambda[1]).
  40 //        Shower is considered as EM if f() > 0 )
  41 // One can visualize current cuts calling method PlotDispersionCuts().
  42 //
  43 // use case:
  44 //  root [0] AliPHOSPIDv0 * p1 = new AliPHOSPIDv0("galice.root")
  45 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  46 //  root [1] p1->SetIdentificationMethod("disp ellipse")
  47 //  root [2] p1->ExecuteTask()
  48 //  root [3] AliPHOSPIDv0 * p2 = new AliPHOSPIDv0("galice1.root","ts1")
  49 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  50 //                // reading headers from file galice1.root and TrackSegments
  51 //                // with title "ts1"
  52 //  root [4] p2->SetRecParticlesBranch("rp1")
  53 //                // set file name for the branch RecParticles
  54 //  root [5] p2->ExecuteTask("deb all time")
  55 //                // available options
  56 //                // "deb" - prints # of reconstructed particles
  57 //                // "deb all" -  prints # and list of RecParticles
  58 //                // "time" - prints benchmarking results
  59 //
  60 //*-- Author: Yves Schutz (SUBATECH)  & Gines Martinez (SUBATECH) &
  61 //            Dmitri Peressounko (SUBATECH & Kurchatov Institute)
  62 //            Completely redesined by Dmitri Peressounko, March 2001
  63
  64 // --- ROOT system ---
  65 #include "TTree.h"
  66 #include "TF2.h"
  67 #include "TFormula.h"
  68 #include "TCanvas.h"
  69
  70 #include "TBenchmark.h"
  71 // --- Standard library ---
  72
  73 // --- AliRoot header files ---
  74 #include "AliLog.h"
  75 #include "AliRun.h"
  76 #include "AliGenerator.h"
  77 #include "AliPHOSPIDv0.h"
  78 #include "AliPHOSEmcRecPoint.h"
  79 #include "AliPHOSTrackSegment.h"
  80 #include "AliPHOSRecParticle.h"
  81 #include "AliPHOSGeometry.h"
  82 #include "AliPHOSLoader.h"
  83
  84 ClassImp( AliPHOSPIDv0)
  85
  86 //____________________________________________________________________________
  87 AliPHOSPIDv0::AliPHOSPIDv0():AliPHOSPID()
  88 {
  89   // default ctor
  90   fFormula           = 0 ;
  91   fDispersion        = 0. ;
  92   fCpvEmcDistance    = 0 ;
  93   fTimeGate          = 2.e-9 ;
  94   fEventFolderName    = "" ;
  95   fTrackSegmentsTitle= "" ;
  96   fRecPointsTitle    = "" ;
  97   fRecParticlesTitle = "" ;
  98   fIDOptions         = "dis time" ;
  99   fRecParticlesInRun = 0 ;
 100   fClusterizer = 0;
 101   fTSMaker = 0;
 102 }
 103
 104 //____________________________________________________________________________
 105 AliPHOSPIDv0::AliPHOSPIDv0(const char * evFolderName,const char * name) : AliPHOSPID(evFolderName, name)
 106 {
 107   //ctor with the indication on where to look for the track segments
 108
 109   fFormula        = new TFormula("LambdaCuts","(x>1)*(x<2.5)*(y>0)*(y<x)") ;
 110   fDispersion     = 2.0 ;
 111   fCpvEmcDistance = 3.0 ;
 112   fTimeGate          = 2.e-9 ;
 113
 114   fEventFolderName     = GetTitle() ;
 115   fTrackSegmentsTitle = GetName();
 116   fRecPointsTitle     = GetName();
 117   fRecParticlesTitle  = GetName();
 118   fIDOptions          = "dis time" ;
 119
 120   fRecParticlesInRun = 0 ;
 121
 122   Init() ;
 123
 124 }
 125
 126 //____________________________________________________________________________
 127 AliPHOSPIDv0::~AliPHOSPIDv0()
 128 {
 129 }
 130
 131 //____________________________________________________________________________
 132 Float_t  AliPHOSPIDv0::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t *  Axis)const
 133 {
 134   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
 135
 136   const AliPHOSGeometry * geom = AliPHOSLoader::GetPHOSGeometry() ;
 137   TVector3 vecEmc ;
 138   TVector3 vecCpv ;
 139
 140   emc->GetLocalPosition(vecEmc) ;
 141   cpv->GetLocalPosition(vecCpv) ;
 142   if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
 143
 144     // Correct to difference in CPV and EMC position due to different distance to center.
 145     // we assume, that particle moves from center
 146     Float_t dCPV = geom->GetIPtoOuterCoverDistance();
 147     Float_t dEMC = geom->GetIPtoCrystalSurface() ;
 148     dEMC         = dEMC / dCPV ;
 149     vecCpv = dEMC * vecCpv  - vecEmc ;
 150     if (Axis == "X") return vecCpv.X();
 151     if (Axis == "Y") return vecCpv.Y();
 152     if (Axis == "Z") return vecCpv.Z();
 153     if (Axis == "R") return vecCpv.Mag();
 154   }
 155
 156   return 100000000 ;
 157 }
 158
 159 //____________________________________________________________________________
 160 void  AliPHOSPIDv0::Exec(Option_t * option)
 161 {
 162   //Steering method
 163
 164   if( strcmp(GetName(), "")== 0 )
 165     Init() ;
 166
 167   if(strstr(option,"tim"))
 168     gBenchmark->Start("PHOSPID");
 169
 170   if(strstr(option,"print")) {
 171     Print() ;
 172     return ;
 173   }
 174
 175   AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
 176   if(runget == 0x0)
 177    {
 178      AliError(Form("Can not find run getter in event folder \"%s\"",
 179                    GetTitle()));
 180      return;
 181    }
 182
 183   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
 184   if ( gime == 0 )
 185    {
 186      AliError("Could not obtain the Loader object !");
 187      return ;
 188    }
 189
 190   if(gime->BranchExists("RecParticles") )
 191     return ;
 192
 193
 194   Int_t nevents = runget->GetNumberOfEvents() ;       //(Int_t) gAlice->TreeE()->GetEntries() ;
 195
 196   Int_t ievent ;
 197
 198   for(ievent = 0; ievent < nevents; ievent++){
 199     runget->GetEvent(ievent);
 200     AliInfo(Form("event %d %d %d",
 201                  ievent, gime->EmcRecPoints(),
 202                  gime->TrackSegments())) ;
 203     MakeRecParticles() ;
 204
 205     WriteRecParticles();
 206
 207     if(strstr(option,"deb"))
 208       PrintRecParticles(option) ;
 209
 210     //increment the total number of rec particles per run
 211     fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ;
 212
 213   }
 214
 215   if(strstr(option,"tim")){
 216     gBenchmark->Stop("PHOSPID");
 217     AliInfo(Form("took %f seconds for PID %f seconds per event",
 218                  gBenchmark->GetCpuTime("PHOSPID"),
 219                  gBenchmark->GetCpuTime("PHOSPID")/nevents)) ;
 220   }
 221
 222 }
 223 //____________________________________________________________________________
 224 void AliPHOSPIDv0::Init()
 225 {
 226   // Make all memory allocations that are not possible in default constructor
 227   // Add the PID task to the list of PHOS tasks
 228
 229   AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
 230   if(runget == 0x0)
 231    {
 232      AliError(Form("Can not find run getter in event folder \"%s\"",
 233                    GetTitle()));
 234      return;
 235    }
 236
 237   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
 238   if ( gime == 0 )
 239    {
 240      AliError("Could not obtain the Loader object !");
 241      return ;
 242    }
 243
 244   gime->PostPID(this);
 245   gime->LoadRecParticles("UPDATE");
 246
 247 }
 248
 249 //____________________________________________________________________________
 250 void  AliPHOSPIDv0::MakeRecParticles()
 251 {
 252   // Reconstructs the particles from the tracksegments
 253
 254   TString taskName(GetName()) ;
 255   taskName.Remove(taskName.Index(Version())-1) ;
 256
 257   AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
 258   if(runget == 0x0)
 259    {
 260      AliError(Form("Can not find run getter in event folder \"%s\"",
 261                    GetTitle()));
 262      return;
 263    }
 264
 265   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
 266   if ( gime == 0 )
 267    {
 268      AliError("Could not obtain the Loader object !");
 269      return ;
 270    }
 271
 272   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 273   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 274   TClonesArray * trackSegments = gime->TrackSegments() ;
 275   TClonesArray * recParticles  = gime->RecParticles() ;
 276
 277   recParticles->Clear();
 278
 279   TIter next(trackSegments) ;
 280   AliPHOSTrackSegment * ts ;
 281   Int_t index = 0 ;
 282   AliPHOSRecParticle * rp ;
 283
 284   Bool_t ellips = fIDOptions.Contains("ell",TString::kIgnoreCase ) ;
 285   Bool_t disp   = fIDOptions.Contains("dis",TString::kIgnoreCase ) ;
 286   Bool_t time   = fIDOptions.Contains("tim",TString::kIgnoreCase ) ;
 287
 288   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
 289
 290     new( (*recParticles)[index] ) AliPHOSRecParticle() ;
 291     rp = (AliPHOSRecParticle *)recParticles->At(index) ;
 292     rp->SetTrackSegment(index) ;
 293     rp->SetIndexInList(index) ;
 294
 295     AliPHOSEmcRecPoint * emc = 0 ;
 296     if(ts->GetEmcIndex()>=0)
 297       emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
 298
 299     AliPHOSRecPoint    * cpv = 0 ;
 300     if(ts->GetCpvIndex()>=0)
 301       cpv = (AliPHOSRecPoint *)   cpvRecPoints->At(ts->GetCpvIndex()) ;
 302
 303     //set momentum and energy first
 304     Float_t    e = emc->GetEnergy() ;
 305     TVector3 dir = GetMomentumDirection(emc,cpv) ;
 306     dir.SetMag(e) ;
 307
 308     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
 309     rp->SetCalcMass(0);
 310
 311     //now set type (reconstructed) of the particle
 312     Int_t showerprofile = 0;  // 0 narrow and 1 wide
 313
 314     if(ellips){
 315       Float_t lambda[2] ;
 316       emc->GetElipsAxis(lambda) ;
 317       if(fFormula->Eval(lambda[0],lambda[1]) <= 0 )
 318         showerprofile = 1 ;  // not narrow
 319     }
 320
 321     if(disp)
 322       if(emc->GetDispersion() > fDispersion )
 323         showerprofile = 1 ;  // not narrow
 324
 325     Int_t slow = 0 ;
 326     if(time)
 327       if(emc->GetTime() > fTimeGate )
 328         slow = 0 ;
 329
 330     // Looking at the CPV detector
 331     Int_t cpvdetector= 0 ;  //1 hit and 0 no hit
 332     if(cpv)
 333       if(GetDistance(emc, cpv,  "R") < fCpvEmcDistance)
 334         cpvdetector = 1 ;
 335
 336     Int_t type = showerprofile + 2 * slow  + 4 * cpvdetector ;
 337     rp->SetType(type) ;
 338     rp->SetProductionVertex(0,0,0,0);
 339     rp->SetFirstMother(-1);
 340     rp->SetLastMother(-1);
 341     rp->SetFirstDaughter(-1);
 342     rp->SetLastDaughter(-1);
 343     rp->SetPolarisation(0,0,0);
 344     index++ ;
 345   }
 346
 347 }
 348
 349 //____________________________________________________________________________
 350 void  AliPHOSPIDv0:: Print(const Option_t *) const
 351 {
 352   // Print the parameters used for the particle type identification
 353   TString message ;
 354   message  = "=============== AliPHOSPIDv0 ================\n" ;
 355   message += "Making PID\n" ;
 356   message += "    Headers file:               %s\n" ;
 357   message += "    RecPoints branch title:     %s\n" ;
 358   message += "    TrackSegments Branch title: %s\n" ;
 359   message += "    RecParticles Branch title   %s\n" ;
 360   message += "with parameters:\n"  ;
 361   message += "    Maximal EMC - CPV  distance (cm) %f\n" ;
 362   AliInfo(Form( message.Data(),
 363        GetTitle(),
 364        fRecPointsTitle.Data(),
 365        fTrackSegmentsTitle.Data(),
 366        fRecParticlesTitle.Data(),
 367        fCpvEmcDistance ));
 368
 369   if(fIDOptions.Contains("dis",TString::kIgnoreCase ))
 370     AliInfo(Form("                    dispersion cut %f",  fDispersion )) ;
 371   if(fIDOptions.Contains("ell",TString::kIgnoreCase ))
 372     AliInfo(Form("             Eliptic cuts function: %s",
 373                  fFormula->GetTitle() )) ;
 374   if(fIDOptions.Contains("tim",TString::kIgnoreCase ))
 375     AliInfo(Form("             Time Gate used: %f",  fTimeGate)) ;
 376 }
 377
 378 //____________________________________________________________________________
 379 void  AliPHOSPIDv0::SetShowerProfileCut(const char * formula)
 380 {
 381   //set shape of the cut on the axis of ellipce, drown around shouer
 382   //shower considered "narrow" if Formula(lambda[0],lambda[1]) > 0.
 383   if(fFormula)
 384     delete fFormula;
 385   fFormula = new TFormula("Lambda Cut",formula) ;
 386 }
 387 //____________________________________________________________________________
 388 void  AliPHOSPIDv0::WriteRecParticles()
 389 {
 390   // Saves the reconstructed particles too a file
 391
 392   AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
 393   if(runget == 0x0)
 394    {
 395      AliError(Form("Can not find run getter in event folder \"%s\"",
 396                    GetTitle()));
 397      return;
 398    }
 399
 400   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
 401   if ( gime == 0 )
 402    {
 403      AliError("Could not obtain the Loader object !");
 404      return ;
 405    }
 406
 407   TClonesArray * recParticles = gime->RecParticles() ;
 408   recParticles->Expand(recParticles->GetEntriesFast() ) ;
 409
 410   TTree * treeR = gime->TreeR();
 411
 412   if(!treeR){
 413     gime->MakeTree("R");
 414     treeR = gime->TreeR() ;
 415   }
 416
 417   //First rp
 418   Int_t bufferSize = 32000 ;
 419   TBranch * rpBranch = treeR->Branch("PHOSRP",&recParticles,bufferSize);
 420   rpBranch->SetTitle(fRecParticlesTitle);
 421
 422   //second, pid
 423   Int_t splitlevel = 0 ;
 424   AliPHOSPIDv0 * pid = this ;
 425   TBranch * pidBranch = treeR->Branch("AliPHOSPID","AliPHOSPIDv0",&pid,bufferSize,splitlevel);
 426   pidBranch->SetTitle(fRecParticlesTitle.Data());
 427
 428   rpBranch->Fill() ;
 429   pidBranch->Fill() ;
 430
 431   gime->WriteRecParticles("OVERWRITE");
 432   gime->WritePID("OVERWRITE");
 433
 434 }
 435
 436 //____________________________________________________________________________
 437 void  AliPHOSPIDv0::PlotDispersionCuts()const
 438 {
 439   // produces a plot of the dispersion cut
 440   TCanvas*  lambdas = new TCanvas("lambdas","Cuts on the ellipse axis",200,10,700,500);
 441
 442   if(fIDOptions.Contains("ell",TString::kIgnoreCase ) ){
 443     TF2 * ell = new TF2("Elliptic Cuts",fFormula->GetName(),0,3,0,3) ;
 444     ell->SetMinimum(0.0000001) ;
 445     ell->SetMaximum(0.001) ;
 446     ell->SetLineStyle(1) ;
 447     ell->SetLineWidth(2) ;
 448     ell->Draw() ;
 449   }
 450
 451   if( fIDOptions.Contains("dis",TString::kIgnoreCase ) ){
 452     TF2 * dsp = new TF2("dispersion","(y<x)*(x*x+y*y < [0]*[0])",0,3,0,3) ;
 453     dsp->SetParameter(0,fDispersion) ;
 454     dsp->SetMinimum(0.0000001) ;
 455     dsp->SetMaximum(0.001) ;
 456     dsp->SetLineStyle(1) ;
 457     dsp->SetLineColor(2) ;
 458     dsp->SetLineWidth(2) ;
 459     dsp->SetNpx(200) ;
 460     dsp->SetNpy(200) ;
 461     if(fIDOptions.Contains("ell",TString::kIgnoreCase ) )
 462       dsp->Draw("same") ;
 463     else
 464       dsp->Draw() ;
 465   }
 466   lambdas->Update();
 467 }
 468
 469 //____________________________________________________________________________
 470 TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const
 471 {
 472   // Calculates the momentum direction:
 473   //   1. if only a EMC RecPoint, direction is given by IP and this RecPoint
 474   //   2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints
 475   //  However because of the poor position resolution of PPSD the direction is always taken as if we were
 476   //  in case 1.
 477
 478   TVector3 dir(0,0,0) ;
 479
 480   TVector3 emcglobalpos ;
 481   TMatrix  dummy ;
 482
 483   emc->GetGlobalPosition(emcglobalpos, dummy) ;
 484
 485
 486   // The following commented code becomes valid once the PPSD provides
 487   // a reasonable position resolution, at least as good as EMC !
 488   //   TVector3 ppsdlglobalpos ;
 489   //   TVector3 ppsduglobalpos ;
 490   //   if( fPpsdLowRecPoint ){ // certainly a photon that has concerted
 491   //     fPpsdLowRecPoint->GetGlobalPosition(ppsdlglobalpos, mdummy) ;
 492   //     dir = emcglobalpos -  ppsdlglobalpos ;
 493   //     if( fPpsdUpRecPoint ){ // not looks like a charged
 494   //        fPpsdUpRecPoint->GetGlobalPosition(ppsduglobalpos, mdummy) ;
 495   //        dir = ( dir +  emcglobalpos -  ppsduglobalpos ) * 0.5 ;
 496   //      }
 497   //   }
 498   //   else { // looks like a neutral
 499   //    dir = emcglobalpos ;
 500   //  }
 501
 502   dir = emcglobalpos ;
 503   dir.SetZ( -dir.Z() ) ;   // why ?
 504   dir.SetMag(1.) ;
 505
 506   //account correction to the position of IP
 507   Float_t xo,yo,zo ; //Coordinates of the origin
 508   gAlice->Generator()->GetOrigin(xo,yo,zo) ;
 509   TVector3 origin(xo,yo,zo);
 510   dir = dir - origin ;
 511
 512   return dir ;
 513 }
 514 //____________________________________________________________________________
 515 void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
 516 {
 517   // Print table of reconstructed particles
 518
 519   AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
 520   if(runget == 0x0)
 521    {
 522      AliError(Form("Can not find run getter in event folder \"%s\"",
 523                    GetTitle()));
 524      return;
 525    }
 526
 527   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
 528   if ( gime == 0 )
 529    {
 530      AliError("Could not obtain the Loader object !");
 531      return ;
 532    }
 533
 534   TString taskName(GetName()) ;
 535   taskName.Remove(taskName.Index(Version())-1) ;
 536   TClonesArray * recParticles = gime->RecParticles() ;
 537
 538   TString message ;
 539   message  = "event %d\n" ;
 540   message += "       found %d RecParticles\n" ;
 541   AliInfo(Form(message.Data(),
 542                gAlice->GetEvNumber(), recParticles->GetEntriesFast() )) ;
 543
 544   if(strstr(option,"all")) {  // printing found TS
 545     AliInfo("  PARTICLE   Index"  ) ;
 546
 547     Int_t index ;
 548     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
 549       AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ;
 550
 551       Text_t particle[11];
 552       switch(rp->GetType()) {
 553       case  AliPHOSFastRecParticle::kNEUTRALEMFAST:
 554         strcpy( particle, "NEUTRAL EM FAST");
 555         break;
 556       case  AliPHOSFastRecParticle::kNEUTRALHAFAST:
 557         strcpy(particle, "NEUTRAL HA FAST");
 558         break;
 559       case  AliPHOSFastRecParticle::kNEUTRALEMSLOW:
 560         strcpy(particle, "NEUTRAL EM SLOW");
 561         break ;
 562       case  AliPHOSFastRecParticle::kNEUTRALHASLOW:
 563         strcpy(particle, "NEUTRAL HA SLOW");
 564         break ;
 565       case  AliPHOSFastRecParticle::kCHARGEDEMFAST:
 566         strcpy(particle, "CHARGED EM FAST") ;
 567         break ;
 568       case  AliPHOSFastRecParticle::kCHARGEDHAFAST:
 569         strcpy(particle, "CHARGED HA FAST") ;
 570         break ;
 571       case  AliPHOSFastRecParticle::kCHARGEDEMSLOW:
 572         strcpy(particle, "CHARGEDEMSLOW") ;
 573         break ;
 574       case  AliPHOSFastRecParticle::kCHARGEDHASLOW:
 575         strcpy(particle, "CHARGED HA SLOW") ;
 576         break ;
 577       }
 578
 579       //    Int_t * primaries;
 580       //    Int_t nprimaries;
 581       //    primaries = rp->GetPrimaries(nprimaries);
 582
 583       AliInfo(Form("          %s     %d",
 584                    particle, rp->GetIndexInList())) ;
 585     }
 586   }
 587 }
 588
 589
 590