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