PHOS/AliPHOSClusterizerv1.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 /* $Log:
  19    1 October 2000. Yuri Kharlov:
  20      AreNeighbours()
  21      PPSD upper layer is considered if number of layers>1
  22
  23    18 October 2000. Yuri Kharlov:
  24      AliPHOSClusterizerv1()
  25      CPV clusterizing parameters added
  26
  27      MakeClusters()
  28      After first PPSD digit remove EMC digits only once
  29 */
  30 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
  31 //////////////////////////////////////////////////////////////////////////////
  32 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and
  33 //  unfolds the clusters having several local maxima.
  34 //  Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
  35 //  PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
  36 //  parameters including input digits branch title, thresholds etc.)
  37 //  This TTask is normally called from Reconstructioner, but can as well be used in
  38 //  standalone mode.
  39 // Use Case:
  40 //  root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname")
  41 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  42 //               // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default")
  43 //               // and saves recpoints in branch named "recpointsname" (default = "digitsname")
  44 //  root [1] cl->ExecuteTask()
  45 //               //finds RecPoints in all events stored in galice.root
  46 //  root [2] cl->SetDigitsBranch("digits2")
  47 //               //sets another title for Digitis (input) branch
  48 //  root [3] cl->SetRecPointsBranch("recp2")
  49 //               //sets another title four output branches
  50 //  root [4] cl->SetEmcLocalMaxCut(0.03)
  51 //               //set clusterization parameters
  52 //  root [5] cl->ExecuteTask("deb all time")
  53 //               //once more finds RecPoints options are
  54 //               // deb - print number of found rec points
  55 //               // deb all - print number of found RecPoints and some their characteristics
  56 //               // time - print benchmarking results
  57
  58 // --- ROOT system ---
  59
  60 #include "TROOT.h"
  61 #include "TFile.h"
  62 #include "TFolder.h"
  63 #include "TMath.h"
  64 #include "TMinuit.h"
  65 #include "TTree.h"
  66 #include "TSystem.h"
  67 #include "TBenchmark.h"
  68
  69 // --- Standard library ---
  70
  71 // --- AliRoot header files ---
  72 #include "AliPHOSClusterizerv1.h"
  73 #include "AliPHOSCpvRecPoint.h"
  74 #include "AliPHOSDigit.h"
  75 #include "AliPHOSDigitizer.h"
  76 #include "AliPHOSEmcRecPoint.h"
  77 #include "AliPHOS.h"
  78 #include "AliPHOSGetter.h"
  79 #include "AliRun.h"
  80 #include "AliPHOSCalibrManager.h"
  81 #include "AliPHOSCalibrationData.h"
  82
  83 ClassImp(AliPHOSClusterizerv1)
  84
  85 //____________________________________________________________________________
  86   AliPHOSClusterizerv1::AliPHOSClusterizerv1() : AliPHOSClusterizer()
  87 {
  88   // default ctor (to be used mainly by Streamer)
  89
  90   InitParameters() ;
  91   fDefaultInit = kTRUE ;
  92 }
  93
  94 //____________________________________________________________________________
  95 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const char* headerFile,const char* name, const Bool_t toSplit)
  96 :AliPHOSClusterizer(headerFile, name, toSplit)
  97 {
  98   // ctor with the indication of the file where header Tree and digits Tree are stored
  99
 100   InitParameters() ;
 101   Init() ;
 102   fDefaultInit = kFALSE ;
 103
 104 }
 105
 106 //____________________________________________________________________________
 107   AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
 108 {
 109   // dtor
 110
 111   delete fPedestals ;
 112   delete fGains ;
 113
 114   fSplitFile = 0 ;
 115
 116 }
 117
 118 //____________________________________________________________________________
 119 const TString AliPHOSClusterizerv1::BranchName() const
 120 {
 121   // Returns branch name of the AliPHOSClusterizer
 122
 123   TString branchName(GetName() ) ;
 124   branchName.Remove(branchName.Index(Version())-1) ;
 125   return branchName ;
 126 }
 127
 128 //____________________________________________________________________________
 129 Float_t  AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId) const
 130 {
 131   // Returns energy value in the cell absId with the digit amplitude amp
 132
 133   if(fPedestals || fGains ){  //use calibration data
 134     if(!fPedestals || !fGains ){
 135       Error("Calibrate","Either Pedestals of Gains not set!") ;
 136       return 0 ;
 137     }
 138     Float_t en=(amp - fPedestals->Data(absId))*fGains->Data(absId) ;
 139     if(en>0)
 140       return en ;
 141     else
 142       return 0. ;
 143   }
 144   else{ //simulation
 145     if(absId <= fEmcCrystals) //calibrate as EMC
 146       return fADCpedestalEmc + amp*fADCchanelEmc ;
 147     else //calibrate as CPV
 148       return fADCpedestalCpv+ amp*fADCchanelCpv ;
 149   }
 150 }
 151
 152 //____________________________________________________________________________
 153 void AliPHOSClusterizerv1::Exec(Option_t * option)
 154 {
 155   // Steering method
 156   if( strcmp(GetName(), "")== 0 )
 157     Init() ;
 158
 159   if(strstr(option,"tim"))
 160     gBenchmark->Start("PHOSClusterizer");
 161
 162   if(strstr(option,"print"))
 163     Print("") ;
 164
 165   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 166   if(gime->BranchExists("RecPoints"))
 167     return ;
 168
 169   //test, if this is real data:
 170   TBranch * eh = gAlice->TreeE()->GetBranch("AliPHOSBeamTestEvent") ;
 171   if(eh){ //Read CalibrationData
 172     AliPHOSCalibrManager * cmngr=AliPHOSCalibrManager::GetInstance() ;
 173     if(!cmngr){ //not defined yet
 174       cmngr=AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
 175     }
 176     if(fPedestals)
 177       delete fPedestals ;
 178     fPedestals = new AliPHOSCalibrationData("Pedestals",fCalibrVersion) ;
 179     cmngr->ReadFromRoot(*fPedestals,fCalibrRun) ;
 180     if(fGains)
 181       delete fGains ;
 182     fGains = new AliPHOSCalibrationData("Gains",fCalibrVersion) ;
 183     cmngr->ReadFromRoot(*fGains,fCalibrRun) ;
 184   }
 185
 186   Int_t nevents = gime->MaxEvent() ;
 187   Int_t ievent ;
 188
 189   for(ievent = 0; ievent < nevents; ievent++){
 190     gime->Event(ievent,"D") ;
 191     Int_t pattern = gime->EventPattern() ;
 192     if(pattern){
 193             printf("pattern %d \n",pattern) ;
 194       if(!fPatterns){
 195         Error("Exec","Patterns for reconstruction of events are not set, use SetPatterns() ") ;
 196         return ;
 197       }
 198       Bool_t skip = kTRUE ;
 199       Int_t npat = fPatterns->GetSize() ;
 200       for(Int_t i = 0; i<npat;i++)
 201         if(pattern==fPatterns->At(i)){
 202           skip = kFALSE ;
 203           break ;
 204         }
 205       if(skip)
 206         continue ;
 207     }
 208
 209     if(ievent == 0)
 210       GetCalibrationParameters() ;
 211
 212     fNumberOfEmcClusters  = fNumberOfCpvClusters  = 0 ;
 213
 214     MakeClusters() ;
 215
 216     if(fToUnfold)
 217       MakeUnfolding() ;
 218
 219     WriteRecPoints(ievent) ;
 220
 221     if(strstr(option,"deb"))
 222       PrintRecPoints(option) ;
 223
 224     //increment the total number of digits per run
 225     fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
 226     fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
 227  }
 228
 229   if(strstr(option,"tim")){
 230     gBenchmark->Stop("PHOSClusterizer");
 231     Info("Exec", "  took %f seconds for Clusterizing %f seconds per event \n",
 232          gBenchmark->GetCpuTime("PHOSClusterizer"),
 233          gBenchmark->GetCpuTime("PHOSClusterizer")/nevents ) ;
 234   }
 235 }
 236
 237 //____________________________________________________________________________
 238 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
 239                                     Int_t nPar, Float_t * fitparameters) const
 240 {
 241   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 242   // The initial values for fitting procedure are set equal to the positions of local maxima.
 243   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 244
 245   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 246   TClonesArray * digits = gime->Digits() ;
 247
 248
 249   gMinuit->mncler();                     // Reset Minuit's list of paramters
 250   gMinuit->SetPrintLevel(-1) ;           // No Printout
 251   gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
 252                                          // To set the address of the minimization function
 253
 254   TList * toMinuit = new TList();
 255   toMinuit->AddAt(emcRP,0) ;
 256   toMinuit->AddAt(digits,1) ;
 257
 258   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 259
 260   // filling initial values for fit parameters
 261   AliPHOSDigit * digit ;
 262
 263   Int_t ierflg  = 0;
 264   Int_t index   = 0 ;
 265   Int_t nDigits = (Int_t) nPar / 3 ;
 266
 267   Int_t iDigit ;
 268
 269   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 270
 271   for(iDigit = 0; iDigit < nDigits; iDigit++){
 272     digit = maxAt[iDigit];
 273
 274     Int_t relid[4] ;
 275     Float_t x = 0.;
 276     Float_t z = 0.;
 277     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 278     geom->RelPosInModule(relid, x, z) ;
 279
 280     Float_t energy = maxAtEnergy[iDigit] ;
 281
 282     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 283     index++ ;
 284     if(ierflg != 0){
 285       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
 286       return kFALSE;
 287     }
 288     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 289     index++ ;
 290     if(ierflg != 0){
 291        Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
 292       return kFALSE;
 293     }
 294     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 295     index++ ;
 296     if(ierflg != 0){
 297       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
 298       return kFALSE;
 299     }
 300   }
 301
 302   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 303                       //  depends on it.
 304   Double_t p1 = 1.0 ;
 305   Double_t p2 = 0.0 ;
 306
 307   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 308   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 309   gMinuit->SetMaxIterations(5);
 310   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 311
 312   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 313
 314   if(ierflg == 4){  // Minimum not found
 315     Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
 316     return kFALSE ;
 317   }
 318   for(index = 0; index < nPar; index++){
 319     Double_t err ;
 320     Double_t val ;
 321     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 322     fitparameters[index] = val ;
 323    }
 324
 325   delete toMinuit ;
 326   return kTRUE;
 327
 328 }
 329
 330 //____________________________________________________________________________
 331 void AliPHOSClusterizerv1::GetCalibrationParameters()
 332 {
 333   // Get calibration parameters from AliPHOSDigitizer
 334
 335   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 336
 337   const TTask * task = gime->Digitizer(BranchName()) ;
 338   if(strcmp(task->IsA()->GetName(),"AliPHOSDigitizer")==0){
 339     const AliPHOSDigitizer * dig = static_cast<const AliPHOSDigitizer *>(task) ;
 340
 341     fADCchanelEmc   = dig->GetEMCchannel() ;
 342     fADCpedestalEmc = dig->GetEMCpedestal();
 343
 344     fADCchanelCpv   = dig->GetCPVchannel() ;
 345     fADCpedestalCpv = dig->GetCPVpedestal() ;
 346   }
 347   else{
 348  ;
 349 //    fCalibrationDB = gime->CalibrationDB();
 350   }
 351 }
 352
 353 //____________________________________________________________________________
 354 void AliPHOSClusterizerv1::Init()
 355 {
 356   // Make all memory allocations which can not be done in default constructor.
 357   // Attach the Clusterizer task to the list of PHOS tasks
 358
 359   if ( strcmp(GetTitle(), "") == 0 )
 360     SetTitle("galice.root") ;
 361
 362   TString branchname = GetName() ;
 363   branchname.Remove(branchname.Index(Version())-1) ;
 364
 365   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(),branchname.Data(), fToSplit ) ;
 366   if ( gime == 0 ) {
 367     Error("Init", "Could not obtain the Getter object !") ;
 368     return ;
 369   }
 370
 371   fSplitFile = 0 ;
 372   if(fToSplit){
 373     // construct the name of the file as /path/EMCAL.SDigits.root
 374     //First - extract full path if necessary
 375     TString fileName(GetTitle()) ;
 376     Ssiz_t islash = fileName.Last('/') ;
 377     if(islash<fileName.Length())
 378       fileName.Remove(islash+1,fileName.Length()) ;
 379     else
 380       fileName="" ;
 381     // Next - append the file name
 382     fileName+="PHOS.RecData." ;
 383     if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
 384       fileName+=branchname ;
 385       fileName+="." ;
 386     }
 387     fileName+="root" ;
 388     // Finally - check if the file already opened or open the file
 389     fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));
 390     if(!fSplitFile)
 391       fSplitFile =  TFile::Open(fileName.Data(),"update") ;
 392   }
 393
 394
 395
 396   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 397   fEmcCrystals = geom->GetNModules() *  geom->GetNCristalsInModule() ;
 398
 399   if(!gMinuit)
 400     gMinuit = new TMinuit(100) ;
 401
 402   gime->PostClusterizer(this) ;
 403   gime->PostRecPoints(branchname) ;
 404
 405 }
 406
 407 //____________________________________________________________________________
 408 void AliPHOSClusterizerv1::InitParameters()
 409 {
 410   // Set initial parameters for clusterization
 411
 412   fNumberOfCpvClusters     = 0 ;
 413   fNumberOfEmcClusters     = 0 ;
 414
 415   fCpvClusteringThreshold  = 0.0;
 416   fEmcClusteringThreshold  = 0.2;
 417
 418   fEmcLocMaxCut            = 0.03 ;
 419   fCpvLocMaxCut            = 0.03 ;
 420
 421   fW0                      = 4.5 ;
 422   fW0CPV                   = 4.0 ;
 423
 424   fEmcTimeGate             = 1.e-8 ;
 425
 426   fToUnfold                = kTRUE ;
 427
 428   fPurifyThreshold         = 0.012 ;
 429
 430   fPedestals = 0 ;
 431   fGains     = 0 ;
 432   fCalibrVersion= "v1" ;
 433   fCalibrRun = 1 ;
 434   fPatterns  = 0 ;
 435
 436   TString clusterizerName( GetName()) ;
 437   if (clusterizerName.IsNull() )
 438     clusterizerName = "Default" ;
 439   clusterizerName.Append(":") ;
 440   clusterizerName.Append(Version()) ;
 441   SetName(clusterizerName) ;
 442   fRecPointsInRun          = 0 ;
 443
 444
 445 }
 446
 447 //____________________________________________________________________________
 448 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
 449 {
 450   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 451   //                                       = 1 are neighbour
 452   //                                       = 2 are not neighbour but do not continue searching
 453   // neighbours are defined as digits having at least a common vertex
 454   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 455   //                                      which is compared to a digit (d2)  not yet in a cluster
 456
 457   const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
 458
 459   Int_t rv = 0 ;
 460
 461   Int_t relid1[4] ;
 462   geom->AbsToRelNumbering(d1->GetId(), relid1) ;
 463
 464   Int_t relid2[4] ;
 465   geom->AbsToRelNumbering(d2->GetId(), relid2) ;
 466
 467   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
 468     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 469     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 470
 471     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
 472       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
 473       rv = 1 ;
 474     }
 475     else {
 476       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
 477         rv = 2; //  Difference in row numbers is too large to look further
 478     }
 479
 480   }
 481   else {
 482
 483     if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
 484       rv=2 ;
 485
 486   }
 487
 488   return rv ;
 489 }
 490
 491
 492 //____________________________________________________________________________
 493 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
 494 {
 495   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
 496
 497   Bool_t rv = kFALSE ;
 498   const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
 499
 500   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
 501
 502   if(digit->GetId() <= nEMC )   rv = kTRUE;
 503
 504   return rv ;
 505 }
 506
 507 //____________________________________________________________________________
 508 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 509 {
 510   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
 511
 512   Bool_t rv = kFALSE ;
 513   const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
 514
 515   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
 516
 517   if(digit->GetId() > nEMC )   rv = kTRUE;
 518
 519   return rv ;
 520 }
 521
 522 //____________________________________________________________________________
 523 void AliPHOSClusterizerv1::WriteRecPoints(Int_t event)
 524 {
 525
 526   // Creates new branches with given title
 527   // fills and writes into TreeR.
 528
 529
 530   AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ;
 531   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 532   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 533   TClonesArray * digits = gime->Digits() ;
 534   TTree * treeR ;
 535
 536   if(fToSplit){
 537     if(!fSplitFile)
 538       return ;
 539     fSplitFile->cd() ;
 540     TString name("TreeR") ;
 541     name += event ;
 542     treeR = dynamic_cast<TTree*>(fSplitFile->Get(name));
 543   }
 544   else{
 545     treeR = gAlice->TreeR();
 546   }
 547
 548   if(!treeR){
 549     if(fSplitFile)
 550       gAlice->MakeTree("R", fSplitFile);
 551     else
 552       gAlice->MakeTree("R", gROOT->GetFile(GetTitle()));
 553     treeR = gAlice->TreeR() ;
 554   }
 555
 556   Int_t index ;
 557   //Evaluate position, dispersion and other RecPoint properties...
 558   for(index = 0; index < emcRecPoints->GetEntries(); index++){
 559     AliPHOSEmcRecPoint * emcrp = static_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(index)) ;
 560     emcrp->Purify(fPurifyThreshold) ;
 561     if(emcrp->GetMultiplicity())
 562       emcrp->EvalAll(fW0,digits) ;
 563     else
 564       emcRecPoints->Remove(emcrp) ;
 565  }
 566
 567   emcRecPoints->Compress() ;
 568   emcRecPoints->Sort() ;
 569   for(index = 0; index < emcRecPoints->GetEntries(); index++)
 570     ((AliPHOSEmcRecPoint *)emcRecPoints->At(index))->SetIndexInList(index) ;
 571
 572   emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
 573
 574   //Now the same for CPV
 575   for(index = 0; index < cpvRecPoints->GetEntries(); index++)
 576     ((AliPHOSRecPoint *)cpvRecPoints->At(index))->EvalAll(fW0CPV,digits)  ;
 577
 578   cpvRecPoints->Sort() ;
 579
 580   for(index = 0; index < cpvRecPoints->GetEntries(); index++)
 581     ((AliPHOSRecPoint *)cpvRecPoints->At(index))->SetIndexInList(index) ;
 582
 583   cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
 584
 585   Int_t bufferSize = 32000 ;
 586   Int_t splitlevel = 0 ;
 587
 588   //First EMC
 589   TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
 590   emcBranch->SetTitle(BranchName());
 591
 592   //Now CPV branch
 593   TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
 594   cpvBranch->SetTitle(BranchName());
 595
 596   //And Finally  clusterizer branch
 597   AliPHOSClusterizerv1 * cl = this ;
 598   TBranch * clusterizerBranch = treeR->Branch("AliPHOSClusterizer","AliPHOSClusterizerv1",
 599                                               &cl,bufferSize,splitlevel);
 600   clusterizerBranch->SetTitle(BranchName());
 601
 602   emcBranch        ->Fill() ;
 603   cpvBranch        ->Fill() ;
 604   clusterizerBranch->Fill() ;
 605
 606   treeR->AutoSave() ;
 607   if(gAlice->TreeR()!=treeR)
 608     treeR->Delete();
 609 }
 610
 611 //____________________________________________________________________________
 612 void AliPHOSClusterizerv1::MakeClusters()
 613 {
 614   // Steering method to construct the clusters stored in a list of Reconstructed Points
 615   // A cluster is defined as a list of neighbour digits
 616
 617   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 618
 619   TObjArray * emcRecPoints = gime->EmcRecPoints(BranchName()) ;
 620   TObjArray * cpvRecPoints = gime->CpvRecPoints(BranchName()) ;
 621   emcRecPoints->Delete() ;
 622   cpvRecPoints->Delete() ;
 623
 624   TClonesArray * digits = gime->Digits() ;
 625    if ( !digits ) {
 626    Fatal("MakeClusters", "Digits with name %s not found !", BranchName().Data() ) ;
 627   }
 628   TClonesArray * digitsC =  (TClonesArray*)digits->Clone() ;
 629
 630
 631   // Clusterization starts
 632
 633   TIter nextdigit(digitsC) ;
 634   AliPHOSDigit * digit ;
 635   Bool_t notremoved = kTRUE ;
 636
 637   while ( (digit = (AliPHOSDigit *)nextdigit()) ) { // scan over the list of digitsC
 638     AliPHOSRecPoint * clu = 0 ;
 639
 640     TArrayI clusterdigitslist(1500) ;
 641     Int_t index ;
 642    // cout << "digit " << digit << " " << Calibrate(digit->GetAmp(),digit->GetId()) << endl ;
 643
 644     if (( IsInEmc (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fEmcClusteringThreshold  ) ||
 645         ( IsInCpv (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fCpvClusteringThreshold  ) ) {
 646       Int_t iDigitInCluster = 0 ;
 647
 648       if  ( IsInEmc(digit) ) {
 649         // start a new EMC RecPoint
 650         if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
 651           emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
 652
 653         emcRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
 654         clu = (AliPHOSEmcRecPoint *) emcRecPoints->At(fNumberOfEmcClusters) ;
 655         fNumberOfEmcClusters++ ;
 656         clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId())) ;
 657         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
 658         iDigitInCluster++ ;
 659         digitsC->Remove(digit) ;
 660
 661       } else {
 662
 663         // start a new CPV cluster
 664         if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
 665           cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
 666
 667         cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
 668
 669         clu =  (AliPHOSCpvRecPoint *) cpvRecPoints->At(fNumberOfCpvClusters)  ;
 670         fNumberOfCpvClusters++ ;
 671         clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId()) ) ;
 672         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;
 673         iDigitInCluster++ ;
 674         digitsC->Remove(digit) ;
 675         nextdigit.Reset() ;
 676
 677         // Here we remove remaining EMC digits, which cannot make a cluster
 678
 679         if( notremoved ) {
 680           while( ( digit = (AliPHOSDigit *)nextdigit() ) ) {
 681             if( IsInEmc(digit) )
 682               digitsC->Remove(digit) ;
 683             else
 684               break ;
 685           }
 686           notremoved = kFALSE ;
 687         }
 688
 689       } // else
 690
 691       nextdigit.Reset() ;
 692
 693       AliPHOSDigit * digitN ;
 694       index = 0 ;
 695       while (index < iDigitInCluster){ // scan over digits already in cluster
 696         digit =  (AliPHOSDigit*)digits->At(clusterdigitslist[index])  ;
 697         index++ ;
 698         while ( (digitN = (AliPHOSDigit *)nextdigit()) ) { // scan over the reduced list of digits
 699           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 700           switch (ineb ) {
 701           case 0 :   // not a neighbour
 702             break ;
 703           case 1 :   // are neighbours
 704             clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->GetId() ) ) ;
 705             clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
 706             iDigitInCluster++ ;
 707             digitsC->Remove(digitN) ;
 708             break ;
 709           case 2 :   // too far from each other
 710             goto endofloop;
 711           } // switch
 712
 713         } // while digitN
 714
 715       endofloop: ;
 716         nextdigit.Reset() ;
 717
 718       } // loop over cluster
 719
 720     } // energy theshold
 721
 722
 723   } // while digit
 724
 725   delete digitsC ;
 726
 727 }
 728
 729 //____________________________________________________________________________
 730 void AliPHOSClusterizerv1::MakeUnfolding()
 731 {
 732   // Unfolds clusters using the shape of an ElectroMagnetic shower
 733   // Performs unfolding of all EMC/CPV clusters
 734
 735   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 736
 737
 738   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 739   TObjArray * emcRecPoints = gime->EmcRecPoints(BranchName()) ;
 740   TObjArray * cpvRecPoints = gime->CpvRecPoints(BranchName()) ;
 741   TClonesArray * digits = gime->Digits() ;
 742
 743   // Unfold first EMC clusters
 744   if(fNumberOfEmcClusters > 0){
 745
 746     Int_t nModulesToUnfold = geom->GetNModules() ;
 747
 748     Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
 749     Int_t index ;
 750     for(index = 0 ; index < numberofNotUnfolded ; index++){
 751
 752       AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint *) emcRecPoints->At(index) ;
 753       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
 754         break ;
 755
 756       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 757       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 758       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 759       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
 760
 761       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 762         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 763         if(emcRecPoint->GetNExMax()!=-1){ //Do not remove clusters which failed to unfold
 764           emcRecPoints->Remove(emcRecPoint);
 765           emcRecPoints->Compress() ;
 766           index-- ;
 767           fNumberOfEmcClusters -- ;
 768           numberofNotUnfolded-- ;
 769         }
 770       }
 771       else
 772         emcRecPoint->SetNExMax(1) ; //Only one local maximum
 773
 774       delete[] maxAt ;
 775       delete[] maxAtEnergy ;
 776     }
 777   }
 778   // Unfolding of EMC clusters finished
 779
 780
 781   // Unfold now CPV clusters
 782   if(fNumberOfCpvClusters > 0){
 783
 784     Int_t nModulesToUnfold = geom->GetNModules() ;
 785
 786     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
 787     Int_t index ;
 788     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
 789
 790       AliPHOSRecPoint * recPoint = (AliPHOSRecPoint *) cpvRecPoints->At(index) ;
 791
 792       if(recPoint->GetPHOSMod()> nModulesToUnfold)
 793         break ;
 794
 795       AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint*) recPoint ;
 796
 797       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 798       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 799       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 800       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
 801
 802       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 803         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 804         cpvRecPoints->Remove(emcRecPoint);
 805         cpvRecPoints->Compress() ;
 806         index-- ;
 807         numberofCpvNotUnfolded-- ;
 808         fNumberOfCpvClusters-- ;
 809       }
 810
 811       delete[] maxAt ;
 812       delete[] maxAtEnergy ;
 813     }
 814   }
 815   //Unfolding of Cpv clusters finished
 816
 817 }
 818
 819 //____________________________________________________________________________
 820 Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t r)
 821 {
 822   // Shape of the shower (see PHOS TDR)
 823   // If you change this function, change also the gradient evaluation in ChiSquare()
 824
 825   Double_t r4    = r*r*r*r ;
 826   Double_t r295  = TMath::Power(r, 2.95) ;
 827   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 828   return shape ;
 829 }
 830
 831 //____________________________________________________________________________
 832 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 833                                                  Int_t nMax,
 834                                                  AliPHOSDigit ** maxAt,
 835                                                  Float_t * maxAtEnergy)
 836 {
 837   // Performs the unfolding of a cluster with nMax overlapping showers
 838
 839   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 840
 841   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 842   const TClonesArray * digits = gime->Digits() ;
 843   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 844   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 845
 846   Int_t nPar = 3 * nMax ;
 847   Float_t * fitparameters = new Float_t[nPar] ;
 848
 849   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
 850   if( !rv ) {
 851     // Fit failed, return
 852     iniEmc->SetNExMax(-1) ;
 853     delete[] fitparameters ;
 854     return ;
 855   }
 856
 857   // create ufolded rec points and fill them with new energy lists
 858   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
 859   // and later correct this number in acordance with actual energy deposition
 860
 861   Int_t nDigits = iniEmc->GetMultiplicity() ;
 862   Float_t * efit = new Float_t[nDigits] ;
 863   Float_t xDigit=0.,zDigit=0.,distance=0. ;
 864   Float_t xpar=0.,zpar=0.,epar=0.  ;
 865   Int_t relid[4] ;
 866   AliPHOSDigit * digit = 0 ;
 867   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 868
 869   Int_t iparam ;
 870   Int_t iDigit ;
 871   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 872     digit = (AliPHOSDigit*) digits->At(emcDigits[iDigit] ) ;
 873     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 874     geom->RelPosInModule(relid, xDigit, zDigit) ;
 875     efit[iDigit] = 0;
 876
 877     iparam = 0 ;
 878     while(iparam < nPar ){
 879       xpar = fitparameters[iparam] ;
 880       zpar = fitparameters[iparam+1] ;
 881       epar = fitparameters[iparam+2] ;
 882       iparam += 3 ;
 883       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 884       distance =  TMath::Sqrt(distance) ;
 885       efit[iDigit] += epar * ShowerShape(distance) ;
 886     }
 887   }
 888
 889
 890   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
 891   // so that energy deposited in each cell is distributed betwin new clusters proportionally
 892   // to its contribution to efit
 893
 894   Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
 895   Float_t ratio ;
 896
 897   iparam = 0 ;
 898   while(iparam < nPar ){
 899     xpar = fitparameters[iparam] ;
 900     zpar = fitparameters[iparam+1] ;
 901     epar = fitparameters[iparam+2] ;
 902     iparam += 3 ;
 903
 904     AliPHOSEmcRecPoint * emcRP = 0 ;
 905
 906     if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
 907
 908       if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
 909         emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
 910
 911       (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
 912       emcRP = (AliPHOSEmcRecPoint *) emcRecPoints->At(fNumberOfEmcClusters);
 913       fNumberOfEmcClusters++ ;
 914       emcRP->SetNExMax((Int_t)nPar/3) ;
 915     }
 916     else{//create new entries in fCpvRecPoints
 917       if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
 918         cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
 919
 920       (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
 921       emcRP = (AliPHOSEmcRecPoint *) cpvRecPoints->At(fNumberOfCpvClusters);
 922       fNumberOfCpvClusters++ ;
 923     }
 924
 925     Float_t eDigit ;
 926     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 927       digit = (AliPHOSDigit*) digits->At( emcDigits[iDigit] ) ;
 928       geom->AbsToRelNumbering(digit->GetId(), relid) ;
 929       geom->RelPosInModule(relid, xDigit, zDigit) ;
 930       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 931       distance =  TMath::Sqrt(distance) ;
 932       ratio = epar * ShowerShape(distance) / efit[iDigit] ;
 933       eDigit = emcEnergies[iDigit] * ratio ;
 934       emcRP->AddDigit( *digit, eDigit ) ;
 935     }
 936   }
 937
 938   delete[] fitparameters ;
 939   delete[] efit ;
 940
 941 }
 942
 943 //_____________________________________________________________________________
 944 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 945 {
 946   // Calculates the Chi square for the cluster unfolding minimization
 947   // Number of parameters, Gradient, Chi squared, parameters, what to do
 948
 949   TList * toMinuit = (TList*) gMinuit->GetObjectFit() ;
 950
 951   AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint*) toMinuit->At(0)  ;
 952   TClonesArray * digits = (TClonesArray*)toMinuit->At(1)  ;
 953
 954
 955
 956   //  AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit
 957
 958   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 959
 960   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
 961
 962   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
 963
 964   const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
 965   fret = 0. ;
 966   Int_t iparam ;
 967
 968   if(iflag == 2)
 969     for(iparam = 0 ; iparam < nPar ; iparam++)
 970       Grad[iparam] = 0 ; // Will evaluate gradient
 971
 972   Double_t efit ;
 973
 974   AliPHOSDigit * digit ;
 975   Int_t iDigit ;
 976
 977   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 978
 979     digit = (AliPHOSDigit*) digits->At( emcDigits[iDigit] ) ;
 980
 981     Int_t relid[4] ;
 982     Float_t xDigit ;
 983     Float_t zDigit ;
 984
 985     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 986
 987     geom->RelPosInModule(relid, xDigit, zDigit) ;
 988
 989      if(iflag == 2){  // calculate gradient
 990        Int_t iParam = 0 ;
 991        efit = 0 ;
 992        while(iParam < nPar ){
 993          Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
 994          iParam++ ;
 995          distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
 996          distance = TMath::Sqrt( distance ) ;
 997          iParam++ ;
 998          efit += x[iParam] * ShowerShape(distance) ;
 999          iParam++ ;
1000        }
1001        Double_t sum=0 ;
1002        if(emcEnergies[iDigit] > 0)
1003         sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ;
1004        // Here we assume, that sigma = sqrt(E)
1005        iParam = 0 ;
1006        while(iParam < nPar ){
1007          Double_t xpar = x[iParam] ;
1008          Double_t zpar = x[iParam+1] ;
1009          Double_t epar = x[iParam+2] ;
1010          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
1011          Double_t shape = sum * ShowerShape(dr) ;
1012          Double_t r4 = dr*dr*dr*dr ;
1013          Double_t r295 = TMath::Power(dr,2.95) ;
1014          Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1015                                          0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1016
1017          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
1018          iParam++ ;
1019          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
1020          iParam++ ;
1021          Grad[iParam] += shape ;                                  // Derivative over energy
1022          iParam++ ;
1023        }
1024      }
1025      efit = 0;
1026      iparam = 0 ;
1027
1028      while(iparam < nPar ){
1029        Double_t xpar = x[iparam] ;
1030        Double_t zpar = x[iparam+1] ;
1031        Double_t epar = x[iparam+2] ;
1032        iparam += 3 ;
1033        Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
1034        distance =  TMath::Sqrt(distance) ;
1035        efit += epar * ShowerShape(distance) ;
1036      }
1037
1038      if(emcEnergies[iDigit] > 0) //for the case if rounding error
1039        fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/
1040                emcEnergies[iDigit] ;
1041      // Here we assume, that sigma = sqrt(E)
1042   }
1043
1044 }
1045
1046 //____________________________________________________________________________
1047 void AliPHOSClusterizerv1::Print(Option_t * option)const
1048 {
1049   // Print clusterizer parameters
1050
1051   TString message ;
1052   TString taskName(GetName()) ;
1053   taskName.ReplaceAll(Version(), "") ;
1054
1055   if( strcmp(GetName(), "") !=0 ) {
1056     // Print parameters
1057     message  = "\n--------------- %s %s -----------\n" ;
1058     message += "Clusterizing digits from the file: %s\n" ;
1059     message += "                           Branch: %s\n" ;
1060     message += "                       EMC Clustering threshold = %f\n" ;
1061     message += "                       EMC Local Maximum cut    = %f\n" ;
1062     message += "                       EMC Logarothmic weight   = %f\n" ;
1063     message += "                       CPV Clustering threshold = %f\n" ;
1064     message += "                       CPV Local Maximum cut    = %f\n" ;
1065     message += "                       CPV Logarothmic weight   = %f\n" ;
1066     if(fToUnfold)
1067       message += " Unfolding on\n" ;
1068     else
1069       message += " Unfolding off\n" ;
1070
1071     message += "------------------------------------------------------------------" ;
1072   }
1073   else
1074     message = " AliPHOSClusterizerv1 not initialized " ;
1075
1076   Info("Print", message.Data(),
1077        taskName.Data(),
1078        GetTitle(),
1079        taskName.Data(),
1080        GetName(),
1081        fEmcClusteringThreshold,
1082        fEmcLocMaxCut,
1083        fW0,
1084        fCpvClusteringThreshold,
1085        fCpvLocMaxCut,
1086        fW0CPV ) ;
1087 }
1088
1089
1090 //____________________________________________________________________________
1091 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1092 {
1093   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1094
1095   TObjArray * emcRecPoints = AliPHOSGetter::GetInstance()->EmcRecPoints(BranchName()) ;
1096   TObjArray * cpvRecPoints = AliPHOSGetter::GetInstance()->CpvRecPoints(BranchName()) ;
1097
1098   if(strstr(option,"deb")){
1099     printf("event %d, found %d EmcRecPoints and %d CPV RecPoints \n",
1100            gAlice->GetEvNumber(),
1101            emcRecPoints->GetEntriesFast(),
1102            cpvRecPoints->GetEntriesFast()) ;
1103   }
1104
1105
1106   if(strstr(option,"all")) {
1107     printf("\nEMC clusters\n" );
1108     printf("Index    Ene(MeV) Multi Module     X     Y   Z  Lambda 1 Lambda 2  # of prim  Primaries list\n" );
1109     Int_t index ;
1110     for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
1111       AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
1112       TVector3  locpos;
1113       rp->GetLocalPosition(locpos);
1114       Float_t lambda[2];
1115       rp->GetElipsAxis(lambda);
1116       Int_t * primaries;
1117       Int_t nprimaries;
1118       primaries = rp->GetPrimaries(nprimaries);
1119       printf("\n%6d  %8.2f  %3d     %2d     %4.1f    %4.1f %4.1f %4f  %4f    %2d     : ",
1120               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1121               locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1122
1123       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1124         printf("%d ", primaries[iprimary] ) ;
1125       }
1126     }
1127     //Now plot CPV recPoints
1128     printf("\n CPV RecPoints \n") ;
1129     printf("Index    Ene(MeV) Module     X     Y   Z  # of prim  Primaries list\n") ;
1130     for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
1131         AliPHOSRecPoint * rp = (AliPHOSRecPoint * )cpvRecPoints->At(index) ;
1132
1133         TVector3  locpos;
1134         rp->GetLocalPosition(locpos);
1135
1136         Int_t * primaries;
1137         Int_t nprimaries ;
1138         primaries = rp->GetPrimaries(nprimaries);
1139         printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f %2d     : ",
1140                 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1141                 locpos.X(), locpos.Y(), locpos.Z(), nprimaries) ;
1142         primaries = rp->GetPrimaries(nprimaries);
1143         for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1144           printf("%d ", primaries[iprimary] ) ;
1145         }
1146       }
1147   }
1148 }
1149