d42b41c162f0775b37620b7a9c20b6ba5b121b0f
[u/mrichter/AliRoot.git] / PHOS / AliPHOSAnalyze.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 // Algorythm class to analyze PHOSv1 events:
20 // Construct histograms and displays them.
21 // Use the macro EditorBar.C for best access to the functionnalities
22 //
23 //*-- Author: Y. Schutz (SUBATECH) & Gines Martinez (SUBATECH)
24 //////////////////////////////////////////////////////////////////////////////
25
26 // --- ROOT system ---
27
28 #include "TFile.h"
29 #include "TH1.h"
30 #include "TPad.h"
31 #include "TH2.h"
32 #include "TParticle.h"
33 #include "TClonesArray.h"
34 #include "TTree.h"
35 #include "TMath.h"
36 #include "TCanvas.h" 
37 #include "TStyle.h" 
38
39 // --- Standard library ---
40
41 #include <iostream.h>
42 #include <stdio.h>
43
44 // --- AliRoot header files ---
45
46 #include "AliRun.h"
47 #include "AliPHOSAnalyze.h"
48 #include "AliPHOSClusterizerv1.h"
49 #include "AliPHOSTrackSegmentMakerv1.h"
50 #include "AliPHOSPIDv1.h"
51 #include "AliPHOSReconstructioner.h"
52 #include "AliPHOSDigit.h"
53 #include "AliPHOSTrackSegment.h"
54 #include "AliPHOSRecParticle.h"
55 #include "AliPHOSIndexToObject.h"
56 #include "AliPHOSCPV.h"
57
58 ClassImp(AliPHOSAnalyze)
59
60
61 //____________________________________________________________________________
62   AliPHOSAnalyze::AliPHOSAnalyze()
63 {
64   // default ctor (useless)
65   
66   fRootFile = 0 ; 
67 }
68
69 //____________________________________________________________________________
70 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
71 {
72   // ctor: analyze events from root file "name"
73   
74   Bool_t ok = OpenRootFile(name)  ; 
75   if ( !ok ) {
76     cout << " AliPHOSAnalyze > Error opening " << name << endl ; 
77   }
78   else { 
79       //========== Get AliRun object from file 
80       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
81
82       //=========== Get the PHOS object and associated geometry from the file      
83       fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
84       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
85  
86       //========== Initializes the Index to Object converter
87       fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ; 
88       //========== Current event number 
89       fEvt = -999 ; 
90
91   }
92   fDebugLevel = 0;
93   ResetHistograms() ;
94 }
95
96 //____________________________________________________________________________
97 AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
98 {
99   // copy ctor
100   ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
101 }
102
103 //____________________________________________________________________________
104 void AliPHOSAnalyze::Copy(TObject & obj)
105 {
106   // copy an analysis into an other one
107   TObject::Copy(obj) ;
108   // I do nothing more because the copy is silly but the Code checkers requires one
109 }
110
111 //____________________________________________________________________________
112 AliPHOSAnalyze::~AliPHOSAnalyze()
113 {
114   // dtor
115
116   if (fRootFile->IsOpen() ) 
117     fRootFile->Close() ; 
118   if(fRootFile) {delete fRootFile ; fRootFile=0 ;}
119   if(fPHOS)     {delete fPHOS     ; fPHOS    =0 ;}
120   if(fClu)      {delete fClu      ; fClu     =0 ;}
121   if(fPID)      {delete fPID      ; fPID     =0 ;}
122   if(fRec)      {delete fRec      ; fRec     =0 ;}
123   if(fTrs)      {delete fTrs      ; fTrs     =0 ;}
124
125 }
126
127 //____________________________________________________________________________
128 void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
129   
130   fhEnergyCorrelations  = new TH2F("hEnergyCorrelations","hEnergyCorrelations",40,  0., 0.15, 30, 0., 3.e-5);
131   //========== Create the Clusterizer
132   fClu = new AliPHOSClusterizerv1() ; 
133   fClu->SetEmcEnergyThreshold(0.05) ; 
134   fClu->SetEmcClusteringThreshold(0.20) ; 
135   fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
136   fClu->SetPpsdClusteringThreshold(0.0000001) ; 
137   fClu->SetLocalMaxCut(0.03) ;
138   fClu->SetCalibrationParameters(0., 0.00000001) ; 
139
140   Int_t ievent;
141   
142   for ( ievent=0; ievent<Nevents; ievent++)
143     {  
144       
145       //========== Event Number>         
146       if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
147         cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;  
148       
149       //=========== Connects the various Tree's for evt
150       gAlice->GetEvent(ievent);
151
152       //=========== Gets the Kine TTree
153       gAlice->TreeK()->GetEvent(0) ;
154       
155       //=========== Get the Digit Tree
156       gAlice->TreeD()->GetEvent(0) ;
157       
158       //========== Creating branches ===================================       
159       AliPHOSRecPoint::RecPointsList ** EmcRecPoints =  fPHOS->EmcRecPoints() ;
160       gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
161       
162       AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
163       gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
164       
165       AliPHOSTrackSegment::TrackSegmentsList **  TrackSegmentsList = fPHOS->TrackSegments() ;
166       if( (*TrackSegmentsList) )
167         (*TrackSegmentsList)->Clear() ;
168       gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
169       
170       AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
171       if( (*RecParticleList) )
172         (*RecParticleList)->Clear() ;
173       gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
174       
175       
176       //=========== Gets the Reconstraction TTree
177       gAlice->TreeR()->GetEvent(0) ;
178             
179       AliPHOSPpsdRecPoint * RecPoint ;
180       Int_t relid[4] ; 
181       TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
182       while( ( RecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) )
183         {
184           if(!(RecPoint->GetUp()) ) {
185             AliPHOSDigit *digit ;
186             Int_t iDigit ;
187             for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++) 
188               {
189                 digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
190                 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;    
191                 if((relid[2]==1)&&(relid[3]==1)&&(relid[0]==RecPoint->GetPHOSMod())){
192                   Float_t ConvertorEnergy = fClu->Calibrate(digit->GetAmp()) ;
193                   fhEnergyCorrelations->Fill(ConvertorEnergy,RecPoint->GetTotalEnergy() );  
194                   break ; 
195                 } 
196               }
197             break ;
198           }
199         }
200     }
201   SaveHistograms() ;
202   fhEnergyCorrelations->Draw("BOX") ;
203 }
204
205 //____________________________________________________________________________
206 void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)    
207 {
208   // analyzes Nevents events in a single PHOS module  
209   // Events should be reconstructed by Reconstruct()
210
211   if ( fRootFile == 0 ) 
212     cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;  
213   else
214     {
215       //========== Booking Histograms
216       cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ; 
217       BookingHistograms();
218
219       Int_t ievent;
220       Int_t relid[4] ; 
221       AliPHOSDigit * digit ;
222       AliPHOSEmcRecPoint * emc ;
223       AliPHOSPpsdRecPoint * ppsd ;
224       //      AliPHOSTrackSegment * tracksegment ;
225       AliPHOSRecParticle * recparticle;
226
227       for ( ievent=0; ievent<Nevents; ievent++)
228         {  
229           //========== Event Number>         
230           if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
231             cout <<  "AnalyzeManyEvents > " << "Event is " << ievent << endl ;  
232
233           //=========== Connects the various Tree's for evt
234           gAlice->GetEvent(ievent);
235
236           //=========== Gets the Digit TTree
237           gAlice->TreeD()->GetEvent(0) ;
238
239           //=========== Gets the number of entries in the Digits array 
240           TIter nextdigit(fPHOS->Digits()) ;
241           while( ( digit = (AliPHOSDigit *)nextdigit() ) )
242             {
243               fGeom->AbsToRelNumbering(digit->GetId(), relid) ;         
244               if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ; 
245               else    
246                 {  
247                   if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp())); 
248                   if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
249                 }
250             }
251           
252
253           //=========== Cluster in module
254           TIter nextEmc(*fPHOS->EmcRecPoints()  ) ;
255           while((emc = (AliPHOSEmcRecPoint *)nextEmc())) 
256             {
257               if ( emc->GetPHOSMod() == module )
258                 {  
259                   fhEmcCluster->Fill(  emc->GetTotalEnergy()  ); 
260                   TIter nextPpsd( *fPHOS->PpsdRecPoints()) ;
261                   while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) 
262                     {
263                       if ( ppsd->GetPHOSMod() == module )
264                         {                         
265                           if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
266                         }
267                     } 
268                 }
269             }
270
271           //=========== Cluster in module PPSD Down
272           TIter nextPpsd(*fPHOS->PpsdRecPoints() ) ;
273           while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) 
274             {
275               if ( ppsd->GetPHOSMod() == module )
276                 { 
277                   if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
278                   if (ppsd->GetUp())  fhVetoCluster     ->Fill(ppsd->GetTotalEnergy()) ;
279                 }
280             }
281
282           //========== TRackSegments in the event
283           TIter nextRecParticle(*fPHOS->RecParticles() ) ; 
284           while((recparticle = (AliPHOSRecParticle *)nextRecParticle())) 
285             {
286               if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
287                 { 
288                   cout << "Particle type is " << recparticle->GetType() << endl ; 
289                   Int_t numberofprimaries = 0 ;
290                   Int_t * listofprimaries = recparticle->GetPrimaries(numberofprimaries) ;
291                   cout << "Number of primaries = " << numberofprimaries << endl ; 
292                   Int_t index ;
293                   for ( index = 0 ; index < numberofprimaries ; index++)
294                     cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;  
295                 }
296             }
297         }   // endfor
298       SaveHistograms();
299     }       // endif
300 }           // endfunction
301
302 //____________________________________________________________________________
303 void AliPHOSAnalyze::ReconstructCPV(Int_t Nevents )    
304 {     
305
306   // Perform reconstruction of EMC and CPV (GPS2 or IHEP) for <Nevents> events
307   // Yuri Kharlov. 19 October 2000
308
309   Int_t ievent ;   
310   for ( ievent=0; ievent<Nevents; ievent++) {  
311     if (ievent==0) {
312       cout << "Analyze > Starting Reconstructing " << endl ; 
313       //========== Create the Clusterizer
314       fClu = new AliPHOSClusterizerv1() ; 
315       fClu->SetEmcEnergyThreshold(0.05) ; 
316       fClu->SetEmcClusteringThreshold(0.20) ; 
317       fClu->SetLocalMaxCut(0.03) ;
318       if      (strcmp(fGeom->GetName(),"GPS2") == 0) {
319         fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
320         fClu->SetPpsdClusteringThreshold(0.0000001) ; 
321       }
322       else if (strcmp(fGeom->GetName(),"IHEP") == 0) {
323         fClu->SetLocalMaxCutCPV(0.03) ;
324         fClu->SetLogWeightCutCPV(4.0) ;
325         fClu->SetPpsdEnergyThreshold    (0.09) ;
326       }
327       fClu->SetCalibrationParameters(0., 0.00000001) ; 
328       
329       //========== Creates the track segment maker
330       fTrs = new AliPHOSTrackSegmentMakerv1()  ;
331       
332       //========== Creates the particle identifier for GPS2 only
333       if      (strcmp(fGeom->GetName(),"GPS2") == 0) {
334         fPID = new AliPHOSPIDv1() ;
335         fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ; 
336       }   
337       
338       //========== Creates the Reconstructioner
339       fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; 
340       if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);     
341     }
342       
343     if (fDebugLevel != 0 ||
344         (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
345       cout <<  "======= Analyze ======> Event " << ievent+1 << endl ;
346     
347     //=========== Connects the various Tree's for evt
348     gAlice->GetEvent(ievent);
349     
350     //=========== Gets the Digit TTree
351     gAlice->TreeD()->GetEvent(0) ;
352     
353     //=========== Do the reconstruction
354     fPHOS->Reconstruction(fRec);
355   }
356
357   if(fClu)      {delete fClu      ; fClu     =0 ;}
358   if(fPID)      {delete fPID      ; fPID     =0 ;}
359   if(fRec)      {delete fRec      ; fRec     =0 ;}
360   if(fTrs)      {delete fTrs      ; fTrs     =0 ;}
361   
362 }
363 //-------------------------------------------------------------------------------------
364
365 void AliPHOSAnalyze::Reconstruct(Int_t Nevents )    
366 {     
367   Int_t ievent ;   
368   for ( ievent=0; ievent<Nevents; ievent++)
369     {  
370       if (ievent==0) 
371         {
372           cout << "Analyze > Starting Reconstructing " << endl ; 
373           //========== Create the Clusterizer
374           fClu = new AliPHOSClusterizerv1() ; 
375           fClu->SetEmcEnergyThreshold(0.05) ; 
376           fClu->SetEmcClusteringThreshold(0.20) ; 
377           fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
378           fClu->SetPpsdClusteringThreshold(0.0000001) ; 
379           fClu->SetLocalMaxCut(0.03) ;
380           fClu->SetCalibrationParameters(0., 0.00000001) ; 
381           
382           //========== Creates the track segment maker
383           fTrs = new AliPHOSTrackSegmentMakerv1()  ;
384           //      fTrs->UnsetUnfoldFlag() ; 
385           
386           //========== Creates the particle identifier
387           fPID = new AliPHOSPIDv1() ;
388           fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ; 
389           
390           //========== Creates the Reconstructioner  
391           fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; 
392 //        fRec -> SetDebugReconstruction(kTRUE);     
393
394         }
395       
396       //========== Event Number>         
397     if ((ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
398       cout <<  "======= Analyze ======> Event " << ievent+1 << endl ;
399       
400       //=========== Connects the various Tree's for evt
401       gAlice->GetEvent(ievent);
402
403       //=========== Gets the Digit TTree
404       gAlice->TreeD()->GetEvent(0) ;
405
406       //=========== Do the reconstruction
407       fPHOS->Reconstruction(fRec);
408     }
409
410   if(fClu)      {delete fClu      ; fClu     =0 ;}
411   if(fPID)      {delete fPID      ; fPID     =0 ;}
412   if(fRec)      {delete fRec      ; fRec     =0 ;}
413   if(fTrs)      {delete fTrs      ; fTrs     =0 ;}
414   
415 }
416 //-------------------------------------------------------------------------------------
417
418 //   TClonesArray AllDigitArray = TClonesArray("AliPHOSDigit",1000) ;
419 //   TClonesArray * PhotonsList ;
420 //   TClonesArray * FalsDigitsList ;
421 //   TClonesArray AllPrimary = TClonesArray("TParticle",5000) ;
422 //   TFile * file2 = new TFile("ph100.root") ; // file with added photons
423 //   gAlice = (AliRun*) file2->Get("gAlice") ;
424 //   Int_t ievent;
425 //   Int_t NDigits[Nevents+1] ;
426 //   NDigits[0]=0 ;
427 //   Int_t NAllDigits = 0;
428 //   Int_t NprimPerEvent = 20 ;
429 //   for (ievent=0; ievent <Nevents; ievent++)
430 //     {
431 //       PhotonsList  = gAlice->Particles();  //Primary
432 //       FalsDigitsList  = ((AliPHOSv1 *)gAlice->GetDetector("PHOS"))->Digits();  //Digits
433 //       gAlice->GetEvent(ievent) ;
434 //       gAlice->TreeD()->GetEvent(0) ;
435 //       gAlice->TreeK()->GetEvent(0) ;
436 //       //Copy Primary
437 //       Int_t Nprim ;
438 //       for(Nprim = 0 ;Nprim < NprimPerEvent ; Nprim++)
439 //      new (AllPrimary[Nprim+ievent*NprimPerEvent])  TParticle(*((TParticle *) PhotonsList->At(Nprim))) ;
440
441 //       //Copy Digits
442 //       TIter nextDigit(FalsDigitsList) ;
443 //       AliPHOSDigit * FalseDigit ;
444 //       NDigits[ievent+1] = NDigits[ievent]+ FalsDigitsList->GetEntriesFast() ; 
445 //       while( (FalseDigit = (AliPHOSDigit *) nextDigit()))
446 //      {        
447 //        new (AllDigitArray[NAllDigits])  AliPHOSDigit(FalseDigit->GetPrimary(1),FalseDigit->GetId(),FalseDigit->GetAmp()) ;
448 //        NAllDigits++ ;
449 //      }         
450 //     }
451 //   file2->Close() ;
452
453
454
455 //       //Add primary particles
456 //       cout << "# of Primaries before add " << PrimaryList->GetEntriesFast() << endl;
457 //      Int_t NTruePrimary = 0 ;  //PrimaryList->GetEntriesFast() ;
458 //       Int_t Nprim ;
459 //       for(Nprim = 0; Nprim < NprimPerEvent; Nprim++)
460 //      new ((*PrimaryList)[NTruePrimary+Nprim])  TParticle(*((TParticle *) AllPrimary.At(Nprim+ievent*NprimPerEvent))) ;
461       
462 //       cout << "# of Primaries after add " << PrimaryList->GetEntriesFast() <<endl;
463
464 //       cout << "Digits before add " << DigitsList->GetEntries() << endl ;
465 //       cout << "Digits to add " <<  NDigits[ievent+1]-  NDigits[ievent]<< endl ;
466       
467       //=========== Add fals digits ==============================
468 //       TIter nextDigit(DigitsList) ;
469 //       AliPHOSDigit * FalseDigit ;
470 //       AliPHOSDigit * RealDigit ;
471 //       Int_t NTrueDigits = DigitsList->GetEntriesFast() ; 
472 //       Int_t Ndigit ;
473 //       for(Ndigit=NDigits[ievent];Ndigit<NDigits[ievent+1];Ndigit++)
474 //      {        
475 //        FalseDigit = (AliPHOSDigit*) AllDigitArray.At(Ndigit) ;
476 //        Bool_t Add = kTRUE ; 
477 //        AliPHOSDigit tmpDigit=AliPHOSDigit(FalseDigit->GetPrimary(1)+NTruePrimary,FalseDigit->GetId(),FalseDigit->GetAmp()) ;
478           
479 //        while( (RealDigit = (AliPHOSDigit *) nextDigit()) && Add)
480 //          {
481 //            if((*RealDigit) == (tmpDigit)) 
482 //              {
483 //                *RealDigit=*RealDigit+tmpDigit ;
484 //                Add = kFALSE ;
485 //              }
486 //          }
487 //        if(Add)
488 //          {
489 //            new ((*DigitsList)[NTrueDigits])  AliPHOSDigit(FalseDigit->GetPrimary(1)+NTruePrimary,FalseDigit->GetId(),FalseDigit->GetAmp()) ;
490 //            ((AliPHOSDigit *)DigitsList->At(NTrueDigits))->SetIndexInList(NTrueDigits) ;
491 //            NTrueDigits++ ;
492 //          }     
493 //      }
494 //       cout << "Digits after add " << DigitsList->GetEntries() << endl ;
495
496
497 //____________________________________________________________________________
498 void AliPHOSAnalyze::ReadAndPrintCPV(Int_t Nevents)
499 {
500   //
501   // Read and print generated and reconstructed hits in CPV
502   // Author: Yuri Kharlov
503   // 12 October 2000
504   //
505
506   cout << "Start CPV Analysis"<< endl ;
507   for ( Int_t ievent=0; ievent<Nevents; ievent++) {  
508       
509     //========== Event Number>         
510     cout << endl <<  "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
511     
512     //=========== Connects the various Tree's for evt
513     gAlice->GetEvent(ievent);
514
515     //=========== Get the Hits Tree
516     gAlice->ResetHits();
517     gAlice->TreeH()->GetEvent(0);
518     
519     //========== Creating branches ===================================
520     AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
521     gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , EmcRecPoints  ) ;
522     
523     AliPHOSRecPoint::RecPointsList ** CpvRecPoints = fPHOS->PpsdRecPoints() ;
524     gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", CpvRecPoints ) ;
525
526     //=========== Gets the Reconstruction TTree
527     gAlice->TreeR()->GetEvent(0) ;
528
529     // Read and print CPV hits
530
531     TClonesArray *CPVhits;
532     for (Int_t iModule=0; iModule < fGeom->GetNModules(); iModule++) {
533       CPVModule   cpvModule = fPHOS->GetCPVModule(iModule);
534       CPVhits   = cpvModule.Hits();
535       Int_t nCPVhits  = CPVhits->GetEntriesFast();
536       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
537         CPVHit        *cpvHit = (CPVHit*)CPVhits->UncheckedAt(ihit);
538         TLorentzVector p      = cpvHit->GetMomentum();
539         Float_t        xgen   = cpvHit->GetX();
540         Float_t        zgen   = cpvHit->GetY();
541         Int_t          ipart  = cpvHit->GetIpart();
542         printf("CPV hit in module %d: ",iModule+1);
543         printf(" p = (%f, %f, %f, %f) GeV,\n",
544                p.Px(),p.Py(),p.Pz(),p.Energy());
545         printf("               xy = (%8.4f, %8.4f) cm, ipart = %d\n",
546                xgen,zgen,ipart);
547       }
548     }
549
550     // Read and print CPV reconstructed points
551
552     TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
553     AliPHOSPpsdRecPoint *cpvRecPoint ;
554     while( ( cpvRecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) ) {
555       TVector3  locpos;
556       cpvRecPoint->GetLocalPosition(locpos);
557       Int_t PHOSModule = cpvRecPoint->GetPHOSMod();
558       printf("CPV recpoint in module %d: (X,Y,Z) = (%f,%f,%f) cm\n",
559              PHOSModule,locpos.X(),locpos.Y(),locpos.Z());
560     }
561   }
562 }
563
564 //____________________________________________________________________________
565 void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
566 {
567   //
568   // Analyzes CPV characteristics
569   // Author: Yuri Kharlov
570   // 9 October 2000
571   //
572
573   // Book histograms
574
575   TH1F *hDx   = new TH1F("hDx"  ,"CPV x-resolution@reconstruction",100,-5. , 5.);
576   TH1F *hDz   = new TH1F("hDz"  ,"CPV z-resolution@reconstruction",100,-5. , 5.);
577   TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",      21,-0.5,20.5);
578
579   cout << "Start CPV Analysis"<< endl ;
580   for ( Int_t ievent=0; ievent<Nevents; ievent++) {  
581       
582     //========== Event Number>         
583     if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
584       cout << endl <<  "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
585     
586     //=========== Connects the various Tree's for evt
587     gAlice->GetEvent(ievent);
588
589     //=========== Get the Hits Tree
590     gAlice->ResetHits();
591     gAlice->TreeH()->GetEvent(0);
592     
593     //========== Creating branches ===================================
594     AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
595     gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , EmcRecPoints  ) ;
596     
597     AliPHOSRecPoint::RecPointsList ** CpvRecPoints = fPHOS->PpsdRecPoints() ;
598     gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", CpvRecPoints ) ;
599
600     //=========== Gets the Reconstruction TTree
601     gAlice->TreeR()->GetEvent(0) ;
602
603     TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
604     AliPHOSEmcRecPoint *cpvRecPoint ;
605     CPVModule           cpvModule;
606     TClonesArray        *CPVhits;
607     while( ( cpvRecPoint = (AliPHOSEmcRecPoint *)nextRP() ) ) {
608       TVector3  locpos;
609       cpvRecPoint->GetLocalPosition(locpos);
610       Int_t PHOSModule = cpvRecPoint->GetPHOSMod();
611       Int_t rpMult     = cpvRecPoint->GetDigitsMultiplicity();
612       Float_t xrec  = locpos.X();
613       Float_t zrec  = locpos.Z();
614       Float_t dxmin = 1.e+10;
615       Float_t dzmin = 1.e+10;
616
617       cpvModule = fPHOS->GetCPVModule(PHOSModule-1);
618       CPVhits   = cpvModule.Hits();
619       Int_t nCPVhits  = CPVhits->GetEntriesFast();
620       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
621         CPVHit        *cpvHit = (CPVHit*)CPVhits->UncheckedAt(ihit);
622         Float_t        xgen   = cpvHit->GetX();
623         Float_t        zgen   = cpvHit->GetY();
624         if ( TMath::Abs(xgen-xrec) < TMath::Abs(dxmin) ) dxmin = xgen-xrec;
625         if ( TMath::Abs(zgen-zrec) < TMath::Abs(dzmin) ) dzmin = zgen-zrec;
626       }
627       cpvModule.Clear();
628       hDx  ->Fill(dxmin);
629       hDz  ->Fill(dzmin);
630       hNrp ->Fill(rpMult);
631     }
632
633   }
634
635   // Save histograms
636
637   Text_t outputname[80] ;
638   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
639   TFile output(outputname,"RECREATE");
640   output.cd();
641
642   hDx  ->Write() ;
643   hDz  ->Write() ;
644   hNrp ->Write() ;
645
646   // Plot histograms
647
648   TCanvas *CPVcanvas = new TCanvas("CPV","CPV analysis",20,20,300,900);
649   gStyle->SetOptStat(111111);
650   gStyle->SetOptFit(1);
651   gStyle->SetOptDate(1);
652   CPVcanvas->Divide(3,1);
653
654   CPVcanvas->cd(1);
655   gPad->SetFillColor(10);
656   hNrp->SetFillColor(16);
657   hNrp->Draw();
658
659   CPVcanvas->cd(2);
660   gPad->SetFillColor(10);
661   hDx->SetFillColor(16);
662   hDx->Fit("gaus");
663   hDx->Draw();
664
665   CPVcanvas->cd(3);
666   gPad->SetFillColor(10);
667   hDz->SetFillColor(16);
668   hDz->Fit("gaus");
669   hDz->Draw();
670
671   CPVcanvas->Print("CPV.ps");
672
673 }
674
675 //____________________________________________________________________________
676  void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )    
677 {
678   // analyzes Nevents events and calculate Energy and Position resolution as well as
679   // probaility of correct indentifiing of the incident particle
680
681   //========== Booking Histograms
682   cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ; 
683   BookResolutionHistograms();
684
685   Int_t Counter[9][5] ;     
686   Int_t i1,i2,TotalInd = 0 ;
687   for(i1 = 0; i1<9; i1++)
688     for(i2 = 0; i2<5; i2++)
689       Counter[i1][i2] = 0 ;
690   
691   Int_t TotalPrimary = 0 ;
692   Int_t TotalRecPart = 0 ;
693   Int_t TotalRPwithPrim = 0 ;
694   Int_t ievent;
695
696   cout << "Start Analysing"<< endl ;
697   for ( ievent=0; ievent<Nevents; ievent++)
698     {  
699       
700       //========== Event Number>         
701       if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
702         cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;  
703       
704       //=========== Connects the various Tree's for evt
705       gAlice->GetEvent(ievent);
706
707       
708
709             //=========== Gets the Kine TTree
710       gAlice->TreeK()->GetEvent(0) ;
711       
712       //=========== Gets the list of Primari Particles
713       TClonesArray * PrimaryList  = gAlice->Particles();     
714
715       TParticle * Primary ;
716       Int_t iPrimary ;
717       for ( iPrimary = 0 ; iPrimary < PrimaryList->GetEntries() ; iPrimary++)
718             {
719               Primary = (TParticle*)PrimaryList->UncheckedAt(iPrimary) ;
720               Int_t PrimaryType = Primary->GetPdgCode() ;
721               if( PrimaryType == 22 ) 
722                 fhPrimary->Fill(Primary->Energy()) ;
723             } 
724       
725       //=========== Get the Digit Tree
726       gAlice->TreeD()->GetEvent(0) ;
727
728       //========== Creating branches ===================================       
729       AliPHOSRecPoint::RecPointsList ** EmcRecPoints =  fPHOS->EmcRecPoints() ;
730       gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
731
732       AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
733       gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
734
735       AliPHOSTrackSegment::TrackSegmentsList **  TrackSegmentsList = fPHOS->TrackSegments() ;
736       if( (*TrackSegmentsList) )
737         (*TrackSegmentsList)->Clear() ;
738       gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
739       
740       AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
741       if( (*RecParticleList) )
742         (*RecParticleList)->Clear() ;
743       gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
744             
745
746       //=========== Gets the Reconstraction TTree
747       gAlice->TreeR()->GetEvent(0) ;
748  
749       cout << ievent << " " << (*EmcRecPoints) << " " <<(*PpsdRecPoints) <<fPHOS->Digits()<< endl ; 
750       cout << "    " << " " << (*EmcRecPoints)->GetEntries() << " " <<(*PpsdRecPoints)->GetEntries() <<fPHOS->Digits()->GetEntries()<< endl ; 
751
752       AliPHOSRecParticle * RecParticle ;
753       Int_t iRecParticle ;
754       for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
755         {
756           RecParticle = (AliPHOSRecParticle *) (*RecParticleList)->At(iRecParticle) ;
757           
758           Int_t ModuleNumberRec ;
759           Double_t RecX, RecZ ;
760           fGeom->ImpactOnEmc(RecParticle->Theta(), RecParticle->Phi(), ModuleNumberRec, RecX, RecZ) ;
761           
762           Double_t MinDistance = 10000 ;
763           Int_t ClosestPrimary = -1 ;
764           
765           Int_t numberofprimaries ;
766           Int_t * listofprimaries  = RecParticle->GetPrimaries(numberofprimaries)  ;
767           Int_t index ;
768           TParticle * Primary ;
769           Double_t Distance = MinDistance ;
770           for ( index = 0 ; index < numberofprimaries ; index++)
771             {
772               Primary = (TParticle*)PrimaryList->UncheckedAt(listofprimaries[index]) ;
773               Int_t ModuleNumber ;
774               Double_t PrimX, PrimZ ;
775               fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ;
776               if(ModuleNumberRec == ModuleNumber)
777                 Distance = TMath::Sqrt((RecX-PrimX)*(RecX-PrimX)+(RecZ-PrimZ)*(RecZ-PrimZ) ) ;
778               if(MinDistance > Distance)
779                 {
780                   MinDistance = Distance ;
781                   ClosestPrimary = listofprimaries[index] ;
782                 }
783             }
784           TotalRecPart++ ;
785           if(ClosestPrimary >=0 )
786             {
787               fhPhotonAllEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ; 
788               fhPhotonAllPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),Distance) ;
789               TotalRPwithPrim++;
790               Int_t PrimaryType = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPdgCode() ;
791               TParticlePDG* PDGparticle = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPDG();
792               Double_t charge =  PDGparticle->Charge() ;
793               Int_t PrimaryCode ;
794               switch(PrimaryType)
795                 {
796                 case 22:
797                   PrimaryCode = 0;  //Photon
798                   break;
799                 case 11 :
800                   PrimaryCode = 1;  //Electron
801                   break;
802                 case -11 :
803                   PrimaryCode = 1;  //positron
804                   break;
805                 case 321 :
806                   PrimaryCode = 4;  //K+
807                   break;
808                 case -321 :
809                   PrimaryCode = 4;  //K-
810                   break;
811                 case 310 :
812                   PrimaryCode = 4;  //K0s
813                   break;
814                 case 130 :
815                   PrimaryCode = 4;  //K0l
816                   break;
817                 default:
818                   if(charge)
819                     PrimaryCode = 2; //Charged hadron
820                   else
821                     PrimaryCode = 3; //Neutral hadron
822                   break;
823                 }
824
825               switch(RecParticle->GetType())
826                 {
827                 case AliPHOSFastRecParticle::kGAMMA:
828                   if(PrimaryType == 22){
829                     fhPhotonEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ; 
830                     fhPhotonPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),Distance) ;
831                     fhPhotonReg->Fill(RecParticle->Energy() ) ;
832                     fhPhotonEM->Fill(RecParticle->Energy() ) ;
833                     fhPhotPhot->Fill(RecParticle->Energy() ) ;
834                   }
835                   if(PrimaryType == 2112){ //neutron
836                     fhNReg->Fill(RecParticle->Energy() ) ;
837                     fhNEM->Fill(RecParticle->Energy() ) ;
838                   }
839                   
840                   if(PrimaryType == -2112){ //neutron ~
841                     fhNBarReg->Fill(RecParticle->Energy() ) ;
842                     fhNBarEM->Fill(RecParticle->Energy() ) ;
843
844                   }
845                   if(PrimaryCode == 2){
846                     fhChargedReg->Fill(RecParticle->Energy() ) ;
847                     fhChargedEM->Fill(RecParticle->Energy() ) ;
848                   }
849
850                   fhAllReg->Fill(RecParticle->Energy() ) ;
851                   fhAllEM->Fill(RecParticle->Energy() ) ;
852                   Counter[0][PrimaryCode]++;
853                   break;
854                 case  AliPHOSFastRecParticle::kELECTRON:
855                   if(PrimaryType == 11 || PrimaryType == -11){
856                     fhElectronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ; 
857                     fhElectronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
858                   }
859                   if(PrimaryType == 22) 
860                     fhPhotElec->Fill(RecParticle->Energy() ) ;
861
862                   Counter[1][PrimaryCode]++;
863                   break;
864                 case  AliPHOSFastRecParticle::kNEUTRALHA:
865                   if(PrimaryType == 22) 
866                     fhPhotNeuH->Fill(RecParticle->Energy() ) ;
867
868                   fhNeutralHadronEnergy->Fill( ((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy()) ; 
869                   fhNeutralHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ,Distance  ) ;
870                   Counter[2][PrimaryCode]++;
871                   break ;
872                 case  AliPHOSFastRecParticle::kNEUTRALEM:
873                   if(PrimaryType == 22 || PrimaryType == 11 || PrimaryType == -11){
874                     fhNeutralEMEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),RecParticle->Energy() ) ; 
875                     fhNeutralEMPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
876                   }
877
878                   if(PrimaryType == 22){ //photon
879                     fhPhotNuEM->Fill(RecParticle->Energy() ) ;
880                     fhPhotonEM->Fill(RecParticle->Energy() ) ;
881                   }
882                   if(PrimaryType == 2112) //neutron
883                     fhNEM->Fill(RecParticle->Energy() ) ;
884                   
885                   if(PrimaryType == -2112) //neutron ~
886                     fhNBarEM->Fill(RecParticle->Energy() ) ;
887
888                   if(PrimaryCode == 2)
889                     fhChargedEM->Fill(RecParticle->Energy() ) ;
890
891                   fhAllEM->Fill(RecParticle->Energy() ) ;
892
893                   Counter[3][PrimaryCode]++;
894                   break ;
895                 case  AliPHOSFastRecParticle::kCHARGEDHA:
896                   if(PrimaryType == 22) //photon
897                     fhPhotChHa->Fill(RecParticle->Energy() ) ;
898                   
899                   fhChargedHadronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),RecParticle->Energy() ) ; 
900                   fhChargedHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
901                   Counter[4][PrimaryCode]++ ;
902                   break ;
903                 case  AliPHOSFastRecParticle::kGAMMAHA:
904                   if(PrimaryType == 22) //photon
905                     fhPhotGaHa->Fill(RecParticle->Energy() ) ;
906                   fhPhotonHadronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy()) ; 
907                   fhPhotonHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
908                   Counter[5][PrimaryCode]++ ;
909                   break ;       
910                 case  AliPHOSFastRecParticle::kABSURDEM:              
911                   Counter[6][PrimaryCode]++ ;
912                   break;
913                 case  AliPHOSFastRecParticle::kABSURDHA:
914                   Counter[7][PrimaryCode]++ ;
915                   break;
916                 default:
917                   Counter[8][PrimaryCode]++ ;
918                   break;
919                 }
920             }
921         }  
922     }   // endfor
923   SaveHistograms();
924   cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
925   cout << "Resolutions: Total primary       " << TotalPrimary << endl ;
926   cout << "Resoluitons: Total reconstracted " << TotalRecPart << endl ;
927   cout << "TotalReconstructed with Primarie " << TotalRPwithPrim << endl ;
928   cout << "                        Primary:   Photon   Electron   Ch. Hadr.  Neutr. Hadr  Kaons" << endl ; 
929   cout << "             Detected as photon       " << Counter[0][0] << "          " << Counter[0][1] << "          " << Counter[0][2] << "          " <<Counter[0][3] << "          " << Counter[0][4] << endl ;
930   cout << "           Detected as electron       " << Counter[1][0] << "          " << Counter[1][1] << "          " << Counter[1][2] << "          " <<Counter[1][3] << "          " << Counter[1][4] << endl ; 
931   cout << "     Detected as neutral hadron       " << Counter[2][0] << "          " << Counter[2][1] << "          " << Counter[2][2] << "          " <<Counter[2][3] << "          " << Counter[2][4] << endl ;
932   cout << "         Detected as neutral EM       " << Counter[3][0] << "          " << Counter[3][1] << "          " << Counter[3][2] << "          " <<Counter[3][3] << "          " << Counter[3][4] << endl ;
933   cout << "     Detected as charged hadron       " << Counter[4][0] << "          " << Counter[4][1] << "          " << Counter[4][2] << "          " <<Counter[4][3] << "          " << Counter[4][4] << endl ;
934   cout << "       Detected as gamma-hadron       " << Counter[5][0] << "          " << Counter[5][1] << "          " << Counter[5][2] << "          " <<Counter[5][3] << "          " << Counter[5][4] << endl ;
935   cout << "          Detected as Absurd EM       " << Counter[6][0] << "          " << Counter[6][1] << "          " << Counter[6][2] << "          " <<Counter[6][3] << "          " << Counter[6][4] << endl ;
936   cout << "      Detected as absurd hadron       " << Counter[7][0] << "          " << Counter[7][1] << "          " << Counter[7][2] << "          " <<Counter[7][3] << "          " << Counter[7][4] << endl ;
937   cout << "          Detected as undefined       " << Counter[8][0] << "          " << Counter[8][1] << "          " << Counter[8][2] << "          " <<Counter[8][3] << "          " << Counter[8][4] << endl ;
938       
939       for(i1 = 0; i1<9; i1++)
940         for(i2 = 0; i2<5; i2++)
941           TotalInd+=Counter[i1][i2] ;
942       cout << "Indentified particles            " << TotalInd << endl ;
943       
944 }           // endfunction
945
946
947 //____________________________________________________________________________
948 void  AliPHOSAnalyze::BookingHistograms()
949 {
950   // Books the histograms where the results of the analysis are stored (to be changed)
951
952   delete fhEmcDigit  ;
953   delete fhVetoDigit  ;
954   delete fhConvertorDigit   ;
955   delete  fhEmcCluster   ;
956   delete fhVetoCluster   ;
957   delete fhConvertorCluster  ;
958   delete fhConvertorEmc  ;
959   
960   fhEmcDigit                = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
961   fhVetoDigit               = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
962   fhConvertorDigit          = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
963   fhEmcCluster              = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
964   fhVetoCluster             = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
965   fhConvertorCluster        = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
966   fhConvertorEmc            = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
967
968 }
969 //____________________________________________________________________________
970 void  AliPHOSAnalyze::BookResolutionHistograms()
971 {
972   // Books the histograms where the results of the Resolution analysis are stored
973
974   if(fhPhotonEnergy)
975     delete fhPhotonEnergy ;
976   if(fhPhotonAllEnergy)
977     delete fhPhotonAllEnergy ;
978   if(fhElectronEnergy)
979     delete fhElectronEnergy ;
980   if(fhElectronAllEnergy)
981     delete fhElectronAllEnergy ;
982   if(fhNeutralHadronEnergy)
983     delete fhNeutralHadronEnergy ;
984   if(fhNeutralEMEnergy)
985     delete fhNeutralEMEnergy ;
986   if(fhNeutralEMAllEnergy)
987     delete fhNeutralEMAllEnergy ;
988   if(fhChargedHadronEnergy)
989     delete fhChargedHadronEnergy ;
990   if(fhPhotonHadronEnergy)
991     delete fhPhotonHadronEnergy ;
992   if(fhPhotonPosition)
993     delete fhPhotonPosition ;
994   if(fhPhotonAllPosition)
995     delete fhPhotonAllPosition ;
996   if(fhElectronPosition)
997     delete fhElectronPosition ;
998   if(fhElectronAllPosition)
999     delete fhElectronAllPosition ;
1000   if(fhNeutralHadronPosition)
1001     delete fhNeutralHadronPosition ;
1002   if(fhNeutralEMPosition)
1003     delete fhNeutralEMPosition ;
1004   if(fhNeutralEMAllPosition)
1005     delete fhNeutralEMAllPosition ;
1006   if(fhChargedHadronPosition)
1007     delete fhChargedHadronPosition ;
1008   if(fhPhotonHadronPosition)
1009     delete fhPhotonHadronPosition ;
1010
1011   fhPhotonEnergy            = new TH2F("hPhotonEnergy",  "hPhotonEnergy",              100, 0., 5., 100, 0., 5.);
1012   fhPhotonAllEnergy         = new TH2F("hPhotonAllEnergy",  "hPhotonAllEnergy",        100, 0., 5., 100, 0., 5.);
1013   fhElectronEnergy          = new TH2F("hElectronEnergy","hElectronEnergy",            100, 0., 5., 100, 0., 5.);
1014   fhElectronAllEnergy       = new TH2F("hElectronAllEnergy","hElectronAllEnergy",      100, 0., 5., 100, 0., 5.);
1015   fhNeutralHadronEnergy     = new TH2F("hNeutralHadronEnergy", "hNeutralHadronEnergy", 100, 0., 5., 100, 0., 5.);
1016   fhNeutralEMEnergy         = new TH2F("hNeutralEMEnergy", "hNeutralEMEnergy",         100, 0., 5., 100, 0., 5.);
1017   fhNeutralEMAllEnergy      = new TH2F("hNeutralEMAllEnergy", "hNeutralEMAllEnergy",   100, 0., 5., 100, 0., 5.);
1018   fhChargedHadronEnergy     = new TH2F("hChargedHadronEnergy", "hChargedHadronEnergy", 100, 0., 5., 100, 0., 5.);
1019   fhPhotonHadronEnergy      = new TH2F("hPhotonHadronEnergy","hPhotonHadronEnergy",    100, 0., 5., 100, 0., 5.);
1020   fhPhotonPosition          = new TH2F("hPhotonPosition","hPhotonPosition",                20, 0., 5., 100, 0., 5.);
1021   fhPhotonAllPosition       = new TH2F("hPhotonAllPosition","hPhotonAllPosition",          20, 0., 5., 100, 0., 5.);
1022   fhElectronPosition        = new TH2F("hElectronPosition","hElectronPosition",            20, 0., 5., 100, 0., 5.);
1023   fhElectronAllPosition     = new TH2F("hElectronAllPosition","hElectronAllPosition",      20, 0., 5., 100, 0., 5.);
1024   fhNeutralHadronPosition   = new TH2F("hNeutralHadronPosition","hNeutralHadronPosition",  20, 0., 5., 100, 0., 5.);
1025   fhNeutralEMPosition       = new TH2F("hNeutralEMPosition","hNeutralEMPosition",          20, 0., 5., 100, 0., 5.);
1026   fhNeutralEMAllPosition    = new TH2F("hNeutralEMAllPosition","hNeutralEMAllPosition",    20, 0., 5., 100, 0., 5.);
1027   fhChargedHadronPosition   = new TH2F("hChargedHadronPosition","hChargedHadronPosition",  20, 0., 5., 100, 0., 5.);
1028   fhPhotonHadronPosition    = new TH2F("hPhotonHadronPosition","hPhotonHadronPosition",    20, 0., 5., 100, 0., 5.);
1029
1030   if(fhPhotonReg)
1031     delete fhPhotonReg ;
1032   if(fhAllReg)
1033     delete fhAllReg ;
1034   if(fhNReg)
1035     delete fhNReg ;
1036   if(fhNReg)
1037     delete fhNReg ;
1038   if(fhNReg)
1039     delete fhNReg ;
1040   
1041   fhPhotonReg = new TH1F("hPhotonReg","hPhotonReg", 20, 0., 5.);
1042   fhAllReg    = new TH1F("hAllReg", "hAllReg",  20, 0., 5.);
1043   fhNReg      = new TH1F("hNReg", "hNReg",  20, 0., 5.);
1044   fhNBarReg   = new TH1F("hNBarReg", "hNBarReg",  20, 0., 5.);
1045   fhChargedReg= new TH1F("hChargedReg", "hChargedReg",  20, 0., 5.);
1046   
1047   if(fhPhotonEM)
1048     delete fhPhotonEM ;
1049   if(fhAllEM)
1050     delete fhAllEM ;
1051   if(fhNEM)
1052     delete fhNEM ;
1053   if(fhNBarEM)
1054     delete fhNBarEM ;
1055   if(fhChargedEM)
1056     delete fhChargedEM ;
1057   
1058   fhPhotonEM = new TH1F("hPhotonEM","hPhotonEM", 20, 0., 5.);
1059   fhAllEM    = new TH1F("hAllEM", "hAllEM",  20, 0., 5.);
1060   fhNEM      = new TH1F("hNEM", "hNEM",  20, 0., 5.);
1061   fhNBarEM   = new TH1F("hNBarEM", "hNBarEM",  20, 0., 5.);
1062   fhChargedEM= new TH1F("hChargedEM", "hChargedEM",  20, 0., 5.);
1063   
1064   if(fhPrimary)
1065     delete fhPrimary ;
1066   fhPrimary= new TH1F("hPrimary", "hPrimary",  20, 0., 5.);
1067
1068   if(fhPhotPhot)
1069     delete fhPhotPhot ;
1070   if(fhPhotElec)
1071     delete fhPhotElec ;
1072   if(fhPhotNeuH)
1073     delete fhPhotNeuH ;
1074   if(fhPhotNuEM)
1075     delete fhPhotNuEM ;
1076   if(fhPhotChHa)
1077     delete fhPhotChHa ;
1078   if(fhPhotGaHa)
1079     delete fhPhotGaHa ;
1080
1081   fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 20, 0., 5.);   //Photon registered as photon
1082   fhPhotElec = new TH1F("hPhotElec","hPhotElec", 20, 0., 5.);   //Photon registered as Electron
1083   fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 20, 0., 5.);   //Photon registered as Neutral Hadron
1084   fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 20, 0., 5.);   //Photon registered as Neutral EM
1085   fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 20, 0., 5.);   //Photon registered as Charged Hadron
1086   fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 20, 0., 5.);   //Photon registered as Gamma-Hadron
1087
1088
1089 }
1090 //____________________________________________________________________________
1091 Bool_t AliPHOSAnalyze::Init(Int_t evt)
1092 {
1093   // Do a few initializations: open the root file
1094   //                           get the AliRun object
1095   //                           defines the clusterizer, tracksegment maker and particle identifier
1096   //                           sets the associated parameters
1097
1098   Bool_t ok = kTRUE ; 
1099   
1100    //========== Open galice root file  
1101
1102   if ( fRootFile == 0 ) {
1103     Text_t * name  = new Text_t[80] ; 
1104     cout << "AnalyzeOneEvent > Enter file root file name : " ;  
1105     cin >> name ; 
1106     Bool_t ok = OpenRootFile(name) ; 
1107     if ( !ok )
1108       cout << " AliPHOSAnalyze > Error opening " << name << endl ; 
1109     else { 
1110       //========== Get AliRun object from file 
1111       
1112       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
1113       
1114       //=========== Get the PHOS object and associated geometry from the file 
1115       
1116       fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
1117       fGeom = fPHOS->GetGeometry();
1118       //      fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
1119
1120     } // else !ok
1121   } // if fRootFile
1122   
1123   if ( ok ) {
1124     
1125     //========== Create the Clusterizer
1126
1127     fClu =  new AliPHOSClusterizerv1() ; 
1128     fClu->SetEmcEnergyThreshold(0.030) ; 
1129     fClu->SetEmcClusteringThreshold(0.20) ; 
1130     fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
1131     fClu->SetPpsdClusteringThreshold(0.0000001) ; 
1132     fClu->SetLocalMaxCut(0.03) ;
1133     fClu->SetCalibrationParameters(0., 0.00000001) ;  
1134     cout <<  "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ; 
1135     fClu->PrintParameters() ; 
1136     
1137     //========== Creates the track segment maker
1138     
1139     fTrs = new AliPHOSTrackSegmentMakerv1() ;
1140     cout <<  "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ; 
1141     //   fTrs->UnsetUnfoldFlag() ;
1142     
1143     //========== Creates the particle identifier
1144     
1145     fPID = new AliPHOSPIDv1() ;
1146     cout <<  "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ; 
1147     //fPID->SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m, Float_t l2M) ; 
1148     fPID->SetShowerProfileCuts(0.7, 2.0 , 0.6 , 1.5) ; 
1149
1150     //========== Creates the Reconstructioner  
1151     
1152     fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
1153 //      fRec -> SetDebugReconstruction(kFALSE);     
1154     fRec -> SetDebugReconstruction(kTRUE);     
1155     
1156     //=========== Connect the various Tree's for evt
1157     
1158     if ( evt == -999 ) {
1159       cout <<  "AnalyzeOneEvent > Enter event number : " ; 
1160       cin >> evt ;  
1161       cout << evt << endl ; 
1162     }
1163     fEvt = evt ; 
1164     gAlice->GetEvent(evt);
1165     
1166     //=========== Get the Digit TTree
1167     
1168     gAlice->TreeD()->GetEvent(0) ;     
1169     
1170   } // ok
1171   
1172   return ok ; 
1173 }
1174
1175
1176 //____________________________________________________________________________
1177 void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
1178 {
1179   // Display particles from the Kine Tree in global Alice (theta, phi) coordinates. 
1180   // One PHOS module at the time.
1181   // The particle type can be selected.
1182   
1183   if (evt == -999) 
1184     evt = fEvt ;
1185
1186   Int_t module ; 
1187   cout <<  "DisplayKineEvent > which module (1-5,  -1: all) ? " ; 
1188   cin >> module ; cout << module << endl ; 
1189
1190   Int_t testparticle ; 
1191   cout << " 22      : PHOTON " << endl 
1192        << " (-)11   : (POSITRON)ELECTRON " << endl 
1193        << " (-)2112 : (ANTI)NEUTRON " << endl  
1194        << " -999    : Everything else " << endl ; 
1195   cout  <<  "DisplayKineEvent > enter PDG particle code to display " ; 
1196   cin >> testparticle ; cout << testparticle << endl ; 
1197
1198   Text_t histoname[80] ;
1199   sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ; 
1200
1201   Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
1202   fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
1203
1204   Double_t theta, phi ; 
1205   fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
1206
1207   Int_t tdim = (Int_t)( (tM - tm) / theta ) ; 
1208   Int_t pdim = (Int_t)( (pM - pm) / phi ) ; 
1209
1210   tm -= theta ; 
1211   tM += theta ; 
1212   pm -= phi ; 
1213   pM += phi ; 
1214
1215   TH2F * histoparticle = new TH2F("histoparticle",  histoname, 
1216                                           pdim, pm, pM, tdim, tm, tM) ; 
1217   histoparticle->SetStats(kFALSE) ;
1218
1219   // Get pointers to Alice Particle TClonesArray
1220
1221   TParticle * particle;
1222   TClonesArray * particlearray  = gAlice->Particles();    
1223
1224   Text_t canvasname[80];
1225   sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
1226   TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ; 
1227
1228   // get the KINE Tree
1229
1230   TTree * kine =  gAlice->TreeK() ; 
1231   Stat_t nParticles =  kine->GetEntries() ; 
1232   cout << "DisplayKineEvent > events in kine " << nParticles << endl ; 
1233   
1234   // loop over particles
1235
1236   Double_t kRADDEG = 180. / TMath::Pi() ; 
1237   Int_t index1 ; 
1238   Int_t nparticlein = 0 ; 
1239   for (index1 = 0 ; index1 < nParticles ; index1++){
1240     Int_t nparticle = particlearray->GetEntriesFast() ;
1241     Int_t index2 ; 
1242     for( index2 = 0 ; index2 < nparticle ; index2++) {         
1243       particle            = (TParticle*)particlearray->UncheckedAt(index2) ;
1244       Int_t  particletype = particle->GetPdgCode() ;
1245       if (testparticle == -999 || testparticle == particletype) { 
1246         Double_t phi        = particle->Phi() ;
1247         Double_t theta      = particle->Theta() ;
1248         Int_t mod ; 
1249         Double_t x, z ; 
1250         fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
1251         if ( mod == module ) {
1252           nparticlein++ ; 
1253           if (particle->Energy() >  fClu->GetEmcClusteringThreshold()  )
1254             histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ; 
1255         } 
1256       } 
1257     }
1258   }
1259   kinecanvas->Draw() ; 
1260   histoparticle->Draw("color") ; 
1261   TPaveText *  pavetext = new TPaveText(294, 100, 300, 101); 
1262   Text_t text[40] ; 
1263   sprintf(text, "Particles: %d ", nparticlein) ;
1264   pavetext->AddText(text) ; 
1265   pavetext->Draw() ; 
1266   kinecanvas->Update(); 
1267
1268 }
1269 //____________________________________________________________________________
1270 void AliPHOSAnalyze::DisplayRecParticles()
1271 {
1272   // Display reconstructed particles in global Alice(theta, phi) coordinates. 
1273   // One PHOS module at the time.
1274   // Click on symbols indicate the reconstructed particle type. 
1275
1276   if (fEvt == -999) {
1277     cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ; 
1278     Text_t answer[1] ; 
1279     cin >> answer ; cout << answer ; 
1280 //     if ( answer == "y" ) 
1281 //       AnalyzeOneEvent() ;
1282   } 
1283     if (fEvt != -999) {
1284       
1285       Int_t module ; 
1286       cout <<  "DisplayRecParticles > which module (1-5,  -1: all) ? " ; 
1287       cin >> module ; cout << module << endl ;
1288       Text_t histoname[80] ; 
1289       sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ; 
1290       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
1291       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
1292       Double_t theta, phi ; 
1293       fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
1294       Int_t tdim = (Int_t)( (tM - tm) / theta ) ; 
1295       Int_t pdim = (Int_t)( (pM - pm) / phi ) ; 
1296       tm -= theta ; 
1297       tM += theta ; 
1298       pm -= phi ; 
1299       TH2F * histoRparticle = new TH2F("histoRparticle",  histoname, 
1300                                        pdim, pm, pM, tdim, tm, tM) ; 
1301       histoRparticle->SetStats(kFALSE) ;
1302       Text_t canvasname[80] ; 
1303       sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
1304       TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ; 
1305       AliPHOSRecParticle::RecParticlesList * rpl = *fPHOS->RecParticles() ; 
1306       Int_t nRecParticles = rpl->GetEntries() ; 
1307       Int_t nRecParticlesInModule = 0 ; 
1308       TIter nextRecPart(rpl) ; 
1309       AliPHOSRecParticle * rp ; 
1310       cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ; 
1311       Double_t kRADDEG = 180. / TMath::Pi() ; 
1312       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
1313         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; 
1314         if ( ts->GetPHOSMod() == module ) {
1315           Int_t numberofprimaries = 0 ;
1316           Int_t * listofprimaries = 0;
1317           rp->GetPrimaries(numberofprimaries) ;
1318           cout << "Number of primaries = " << numberofprimaries << endl ; 
1319           Int_t index ;
1320           for ( index = 0 ; index < numberofprimaries ; index++)
1321             cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;  
1322           
1323           nRecParticlesInModule++ ; 
1324           Double_t theta = rp->Theta() * kRADDEG ;
1325           Double_t phi   = rp->Phi() * kRADDEG ;
1326           Double_t energy = rp->Energy() ; 
1327           histoRparticle->Fill(phi, theta, energy) ;
1328         }
1329       }
1330       histoRparticle->Draw("color") ; 
1331
1332       nextRecPart.Reset() ; 
1333       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
1334         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; 
1335         if ( ts->GetPHOSMod() == module )  
1336           rp->Draw("P") ; 
1337       }
1338
1339       Text_t text[80] ; 
1340       sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
1341       TPaveText *  pavetext = new TPaveText(292, 100, 300, 101); 
1342       pavetext->AddText(text) ; 
1343       pavetext->Draw() ; 
1344       rparticlecanvas->Update() ; 
1345     }
1346 }
1347
1348 //____________________________________________________________________________
1349 void AliPHOSAnalyze::DisplayRecPoints()
1350 {
1351   // Display reconstructed points in local PHOS-module (x, z) coordinates. 
1352   // One PHOS module at the time.
1353   // Click on symbols displays the EMC cluster, or PPSD information.
1354
1355   if (fEvt == -999) {
1356     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ; 
1357     Text_t answer[1] ; 
1358     cin >> answer ; cout << answer ; 
1359 //     if ( answer == "y" ) 
1360 //       AnalyzeOneEvent() ;
1361   } 
1362     if (fEvt != -999) {
1363       
1364       Int_t module ; 
1365       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ; 
1366       cin >> module ; cout << module << endl ; 
1367
1368       Text_t canvasname[80];
1369       sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
1370       TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ; 
1371       modulecanvas->Draw() ;
1372
1373       //=========== Creating 2d-histogram of the PHOS module
1374       // a little bit junkie but is used to test Geom functinalities
1375
1376       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
1377       
1378       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); 
1379       // convert angles into coordinates local to the EMC module of interest
1380
1381       Int_t emcModuleNumber ;
1382       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
1383       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
1384       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
1385       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
1386       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;  
1387       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
1388       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; 
1389       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; 
1390       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; 
1391       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;     
1392       Text_t histoname[80];
1393       sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
1394       TH2F * hModule = new TH2F("HistoReconstructed", histoname,
1395                                 xdim, xmin, xMax, zdim, zmin, zMax) ;  
1396       hModule->SetMaximum(2.0);
1397       hModule->SetMinimum(0.0);
1398       hModule->SetStats(kFALSE); 
1399
1400       TIter next(fPHOS->Digits()) ;
1401       Float_t energy, y, z;
1402       Float_t etot=0.;
1403       Int_t relid[4]; Int_t nDigits = 0 ;
1404       AliPHOSDigit * digit ; 
1405
1406       // Making 2D histogram of the EMC module
1407       while((digit = (AliPHOSDigit *)next())) 
1408         {  
1409           fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
1410           if (relid[0] == module && relid[1] == 0)  
1411             {  
1412               energy = fClu->Calibrate(digit->GetAmp()) ;
1413               cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl; 
1414               if (energy >  fClu->GetEmcEnergyThreshold()  ){
1415                 nDigits++ ;
1416                 etot += energy ; 
1417                 fGeom->RelPosInModule(relid,y,z) ;   
1418                 hModule->Fill(y, z, energy) ;
1419               }
1420             } 
1421         }
1422       cout <<"DrawRecPoints >  Found in module " 
1423            << module << " " << nDigits << "  digits with total energy " << etot << endl ;
1424       hModule->Draw("col2") ;
1425
1426       //=========== Cluster in module
1427
1428       //      TClonesArray * emcRP = fPHOS->EmcClusters() ; 
1429       TObjArray * emcRP = *(fPHOS->EmcRecPoints()) ; 
1430       
1431       etot = 0.; 
1432       Int_t totalnClusters = 0 ; 
1433       Int_t nClusters = 0 ;
1434       TIter nextemc(emcRP) ;
1435       AliPHOSEmcRecPoint * emc ;
1436       while((emc = (AliPHOSEmcRecPoint *)nextemc())) 
1437         {
1438           //      Int_t numberofprimaries ;
1439           //      Int_t * primariesarray = new Int_t[10] ;
1440           //      emc->GetPrimaries(numberofprimaries, primariesarray) ;
1441           totalnClusters++ ;
1442           if ( emc->GetPHOSMod() == module )
1443             { 
1444               nClusters++ ; 
1445               energy = emc->GetTotalEnergy() ;   
1446               etot+= energy ;  
1447               emc->Draw("M") ;
1448             }
1449         }
1450       cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ; 
1451       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " EMC Clusters " << endl ;
1452       cout << "DrawRecPoints > total energy  " << etot << endl ; 
1453
1454       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90); 
1455       Text_t text[40] ; 
1456       sprintf(text, "digits: %d;  clusters: %d", nDigits, nClusters) ;
1457       pavetext->AddText(text) ; 
1458       pavetext->Draw() ; 
1459       modulecanvas->Update(); 
1460  
1461       //=========== Cluster in module PPSD Down
1462
1463       //      TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
1464       TObjArray * ppsdRP = *(fPHOS->PpsdRecPoints() );
1465  
1466       etot = 0.; 
1467       TIter nextPpsd(ppsdRP) ;
1468       AliPHOSPpsdRecPoint * ppsd ;
1469       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) 
1470         {
1471           totalnClusters++ ;
1472           if ( ppsd->GetPHOSMod() == module )
1473             { 
1474               nClusters++ ; 
1475               energy = ppsd->GetEnergy() ;   
1476               etot+=energy ;  
1477               if (!ppsd->GetUp()) ppsd->Draw("P") ;
1478             }
1479         }
1480       cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ; 
1481       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Down Clusters " << endl ;
1482       cout << "DrawRecPoints > total energy  " << etot << endl ; 
1483
1484       //=========== Cluster in module PPSD Up
1485   
1486       ppsdRP = *(fPHOS->PpsdRecPoints()) ;
1487      
1488       etot = 0.; 
1489       TIter nextPpsdUp(ppsdRP) ;
1490       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp())) 
1491         {
1492           totalnClusters++ ;
1493           if ( ppsd->GetPHOSMod() == module )
1494             { 
1495               nClusters++ ; 
1496               energy = ppsd->GetEnergy() ;   
1497               etot+=energy ;  
1498               if (ppsd->GetUp()) ppsd->Draw("P") ;
1499             }
1500         }
1501   cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ; 
1502   cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Up Clusters " << endl ;
1503   cout << "DrawRecPoints > total energy  " << etot << endl ; 
1504     
1505     } // if !-999
1506 }
1507
1508 //____________________________________________________________________________
1509 void AliPHOSAnalyze::DisplayTrackSegments()
1510 {
1511   // Display track segments in local PHOS-module (x, z) coordinates. 
1512   // One PHOS module at the time.
1513   // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD.
1514
1515   if (fEvt == -999) {
1516     cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ; 
1517     Text_t answer[1] ; 
1518     cin >> answer ; cout << answer ; 
1519 //     if ( answer == "y" ) 
1520 //       AnalyzeOneEvent() ;
1521   } 
1522     if (fEvt != -999) {
1523
1524       Int_t module ; 
1525       cout <<  "DisplayTrackSegments > which module (1-5,  -1: all) ? " ; 
1526       cin >> module ; cout << module << endl ; 
1527       //=========== Creating 2d-histogram of the PHOS module
1528       // a little bit junkie but is used to test Geom functinalities
1529       
1530       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
1531       
1532       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); 
1533       // convert angles into coordinates local to the EMC module of interest
1534       
1535       Int_t emcModuleNumber ;
1536       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
1537       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
1538       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
1539       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
1540       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;  
1541       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
1542       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; 
1543       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; 
1544       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; 
1545       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;     
1546       Text_t histoname[80];
1547       sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ; 
1548       TH2F * histotrack = new TH2F("histotrack",  histoname, 
1549                                    xdim, xmin, xMax, zdim, zmin, zMax) ;  
1550       histotrack->SetStats(kFALSE); 
1551       Text_t canvasname[80];
1552       sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
1553       TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ; 
1554       histotrack->Draw() ; 
1555
1556       AliPHOSTrackSegment::TrackSegmentsList * trsegl = *(fPHOS->TrackSegments()) ;
1557       AliPHOSTrackSegment * trseg ;
1558  
1559       Int_t nTrackSegments = trsegl->GetEntries() ;
1560       Int_t index ;
1561       Float_t etot = 0 ;
1562       Int_t nTrackSegmentsInModule = 0 ; 
1563       for(index = 0; index < nTrackSegments ; index++){
1564         trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
1565         etot+= trseg->GetEnergy() ;
1566         if ( trseg->GetPHOSMod() == module ) { 
1567           nTrackSegmentsInModule++ ; 
1568           trseg->Draw("P");
1569         }
1570       } 
1571       Text_t text[80] ; 
1572       sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
1573       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90); 
1574       pavetext->AddText(text) ; 
1575       pavetext->Draw() ; 
1576       trackcanvas->Update() ; 
1577       cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
1578     
1579    }
1580 }
1581 //____________________________________________________________________________
1582 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1583 {
1584   // Open the root file named "name"
1585   
1586   fRootFile   = new TFile(name, "update") ;
1587   return  fRootFile->IsOpen() ; 
1588 }
1589 //____________________________________________________________________________
1590 // void AliPHOSAnalyze::SavingHistograms()
1591 // {
1592 //   // Saves the histograms in a root file named "name.analyzed" 
1593
1594 //   Text_t outputname[80] ;
1595 //   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1596 //   TFile output(outputname,"RECREATE");
1597 //   output.cd();
1598 //   if (fhEmcDigit )         
1599 //     fhEmcDigit->Write()  ;
1600 //   if (fhVetoDigit )  
1601 //     fhVetoDigit->Write()  ;
1602 //   if (fhConvertorDigit ) 
1603 //     fhConvertorDigit->Write()   ;
1604 //   if (fhEmcCluster   )
1605 //     fhEmcCluster->Write()   ;
1606 //   if (fhVetoCluster ) 
1607 //     fhVetoCluster->Write()   ;
1608 //   if (fhConvertorCluster )
1609 //     fhConvertorCluster->Write()  ;
1610 //   if (fhConvertorEmc ) 
1611 //     fhConvertorEmc->Write()  ;
1612 //   if (fhPhotonEnergy)    
1613 //     fhPhotonEnergy->Write() ;
1614 //   if (fhPhotonPositionX)  
1615 //     fhPhotonPositionX->Write() ;
1616 //   if (fhPhotonPositionY)  
1617 //     fhPhotonPositionX->Write() ;
1618 //   if (fhElectronEnergy)  
1619 //     fhElectronEnergy->Write() ;
1620 //   if (fhElectronPositionX)
1621 //     fhElectronPositionX->Write() ;
1622 //   if (fhElectronPositionY) 
1623 //     fhElectronPositionX->Write() ;
1624 //   if (fhNeutralHadronEnergy) 
1625 //     fhNeutralHadronEnergy->Write() ;
1626 //   if (fhNeutralHadronPositionX)
1627 //     fhNeutralHadronPositionX->Write() ;
1628 //   if (fhNeutralHadronPositionY) 
1629 //     fhNeutralHadronPositionX->Write() ;
1630 //   if (fhNeutralEMEnergy)   
1631 //     fhNeutralEMEnergy->Write() ;
1632 //   if (fhNeutralEMPositionX)
1633 //     fhNeutralEMPositionX->Write() ;
1634 //   if (fhNeutralEMPositionY) 
1635 //     fhNeutralEMPositionX->Write() ;
1636 //   if (fhChargedHadronEnergy) 
1637 //     fhChargedHadronEnergy->Write() ;
1638 //   if (fhChargedHadronPositionX) 
1639 //     fhChargedHadronPositionX->Write() ;
1640 //   if (fhChargedHadronPositionY)
1641 //     fhChargedHadronPositionX->Write() ;
1642 //   if (fhPhotonHadronEnergy) 
1643 //     fhPhotonHadronEnergy->Write() ;
1644 //   if (fhPhotonHadronPositionX) 
1645 //     fhPhotonHadronPositionX->Write() ;
1646 //   if (fhPhotonHadronPositionY)
1647 //     fhPhotonHadronPositionX->Write() ;
1648
1649 //   output.Write();
1650 //   output.Close();
1651 // }
1652 //____________________________________________________________________________
1653 void AliPHOSAnalyze::SaveHistograms()
1654 {
1655   // Saves the histograms in a root file named "name.analyzed" 
1656
1657   Text_t outputname[80] ;
1658   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1659   TFile output(outputname,"RECREATE");
1660   output.cd();
1661
1662   if (fhPhotonEnergy)    
1663     fhPhotonEnergy->Write() ;
1664   if (fhPhotonAllEnergy)    
1665     fhPhotonAllEnergy->Write() ;
1666   if (fhPhotonPosition)  
1667     fhPhotonPosition->Write() ;
1668   if (fhPhotonAllPosition)  
1669     fhPhotonAllPosition->Write() ;
1670   if (fhElectronEnergy)  
1671     fhElectronEnergy->Write() ;
1672   if (fhElectronAllEnergy)  
1673     fhElectronAllEnergy->Write() ;
1674   if (fhElectronPosition)
1675     fhElectronPosition->Write() ;
1676   if (fhElectronAllPosition)
1677     fhElectronAllPosition->Write() ;
1678   if (fhNeutralHadronEnergy) 
1679     fhNeutralHadronEnergy->Write() ;
1680   if (fhNeutralHadronPosition)
1681     fhNeutralHadronPosition->Write() ;
1682   if (fhNeutralEMEnergy)   
1683     fhNeutralEMEnergy->Write() ;
1684   if (fhNeutralEMAllEnergy)   
1685     fhNeutralEMAllEnergy->Write() ;
1686   if (fhNeutralEMPosition)
1687     fhNeutralEMPosition->Write() ;
1688   if (fhNeutralEMAllPosition)
1689     fhNeutralEMAllPosition->Write() ;
1690   if (fhChargedHadronEnergy) 
1691     fhChargedHadronEnergy->Write() ;
1692   if (fhChargedHadronPosition) 
1693     fhChargedHadronPosition->Write() ;
1694   if (fhPhotonHadronEnergy) 
1695     fhPhotonHadronEnergy->Write() ;
1696   if (fhPhotonHadronPosition) 
1697     fhPhotonHadronPosition->Write() ;
1698   if (fhPhotonReg) 
1699     fhPhotonReg->Write() ;
1700   if (fhAllReg) 
1701     fhAllReg->Write() ;
1702   if(fhNReg)
1703     fhNReg->Write() ;
1704   if(fhNBarReg)
1705     fhNBarReg->Write() ;
1706   if(fhChargedReg)
1707     fhChargedReg->Write() ;
1708   if (fhPhotonEM) 
1709     fhPhotonEM->Write() ;
1710   if (fhAllEM) 
1711     fhAllEM->Write() ;
1712   if(fhNEM)
1713     fhNEM->Write() ;
1714   if(fhNBarEM)
1715     fhNBarEM->Write() ;
1716   if(fhChargedEM)
1717     fhChargedEM->Write() ;
1718   if(fhPrimary)
1719     fhPrimary->Write() ;
1720   if(fhPhotPhot)
1721     fhPhotPhot->Write() ;
1722   if(fhPhotElec)
1723     fhPhotElec->Write() ;
1724   if(fhPhotNeuH)
1725     fhPhotNeuH->Write() ;
1726   if(fhPhotNuEM)
1727     fhPhotNuEM->Write() ;
1728   if(fhPhotNuEM)
1729     fhPhotNuEM->Write() ;
1730   if(fhPhotChHa)
1731     fhPhotChHa->Write() ;
1732   if(fhPhotGaHa)
1733     fhPhotGaHa->Write() ;
1734   if(fhEnergyCorrelations)
1735     fhEnergyCorrelations->Write() ;
1736   
1737   output.Write();
1738   output.Close();
1739 }
1740 //____________________________________________________________________________
1741 void AliPHOSAnalyze::ResetHistograms()
1742 {
1743    fhEnergyCorrelations = 0 ;     //Energy correlations between Eloss in Convertor and PPSD(2)
1744
1745    fhEmcDigit = 0 ;               // Histo of digit energies in the Emc 
1746    fhVetoDigit = 0 ;              // Histo of digit energies in the Veto 
1747    fhConvertorDigit = 0 ;         // Histo of digit energies in the Convertor
1748    fhEmcCluster = 0 ;             // Histo of Cluster energies in Emc
1749    fhVetoCluster = 0 ;            // Histo of Cluster energies in Veto
1750    fhConvertorCluster = 0 ;       // Histo of Cluster energies in Convertor
1751    fhConvertorEmc = 0 ;           // 2d Convertor versus Emc energies
1752
1753    fhPhotonEnergy = 0 ;           // Spectrum of detected photons with photon primary
1754    fhPhotonAllEnergy = 0 ;        // Total spectrum of detected photons
1755    fhElectronEnergy = 0 ;         // Spectrum of detected electrons with electron primary
1756    fhElectronAllEnergy = 0 ;      // Total spectrum of detected electrons
1757    fhNeutralHadronEnergy = 0 ;    // Spectrum of detected neutral hadron
1758    fhNeutralEMEnergy = 0 ;        // Spectrum of detected neutral EM with EM primary
1759    fhNeutralEMAllEnergy = 0 ;     // Spectrum of detected neutral EM
1760    fhChargedHadronEnergy = 0 ;    // Spectrum of detected charged
1761    fhPhotonHadronEnergy = 0 ;     // Spectrum of detected Photon-Hadron
1762    fhPhotonPosition = 0 ;        // Position Resolution of  photons with photon primary
1763    fhPhotonAllPosition = 0 ;     // Position Resolution of  photons
1764    fhElectronPosition = 0 ;      // Position Resolution of electrons with electron primary
1765    fhElectronAllPosition = 0 ;   // Position Resolution of electrons
1766    fhNeutralHadronPosition = 0 ; // Position Resolution of neutral hadron
1767    fhNeutralEMPosition = 0 ;     // Position Resolution of neutral EM with EM primary
1768    fhNeutralEMAllPosition = 0 ;  // Position Resolution of neutral EM
1769    fhChargedHadronPosition = 0 ; // Position Resolution of charged
1770    fhPhotonHadronPosition = 0 ;  // Position Resolution of Photon-Hadron
1771    fhPhotonPositionY = 0 ;        // Y distribution of detected photons
1772    fhElectronPositionY = 0 ;      // Y distribution of detected electrons
1773    fhNeutralHadronPositionY = 0 ; // Y distribution of detected neutral hadron
1774    fhNeutralEMPositionY = 0 ;     // Y distribution of detected neutral EM
1775    fhChargedHadronPositionY = 0 ; // Y distribution of detected charged
1776    fhPhotonHadronPositionY = 0 ;  // Y distribution of detected Photon-Hadron
1777    fhPhotonReg = 0 ;          
1778    fhAllReg = 0 ;          
1779    fhNReg = 0 ;          
1780    fhNBarReg = 0 ;          
1781    fhChargedReg = 0 ;          
1782    fhPhotonEM = 0 ;          
1783    fhAllEM = 0 ;          
1784    fhNEM = 0 ;          
1785    fhNBarEM = 0 ;          
1786    fhChargedEM = 0 ;          
1787    fhPrimary = 0 ;          
1788
1789    fhPhotPhot = 0 ;
1790    fhPhotElec = 0 ;
1791    fhPhotNeuH = 0 ;
1792    fhPhotNuEM = 0 ; 
1793    fhPhotChHa = 0 ;
1794    fhPhotGaHa = 0 ;
1795
1796
1797 }
1798
1799