5c1acbd5ceb0b7beda3676543996c3f0db8b7e28
[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 "AliPHOSHit.h"
57 #include "AliPHOSCPVHit.h"
58 #include "AliPHOSCpvRecPoint.h"
59
60 ClassImp(AliPHOSAnalyze)
61
62 //____________________________________________________________________________
63   AliPHOSAnalyze::AliPHOSAnalyze()
64 {
65   // default ctor (useless)
66   
67   fRootFile = 0 ; 
68 }
69
70 //____________________________________________________________________________
71 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
72 {
73   // ctor: analyze events from root file "name"
74   
75   Bool_t ok = OpenRootFile(name)  ; 
76   if ( !ok ) {
77     cout << " AliPHOSAnalyze > Error opening " << name << endl ; 
78   }
79   else { 
80       //========== Get AliRun object from file 
81       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
82
83       //=========== Get the PHOS object and associated geometry from the file      
84       fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
85       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
86  
87       //========== Initializes the Index to Object converter
88       fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ; 
89       //========== Current event number 
90       fEvt = -999 ; 
91
92   }
93   fDebugLevel = 0;
94   fClu = 0 ;
95   fPID = 0 ;
96   fTrs = 0 ;
97   fRec = 0 ;
98   ResetHistograms() ;
99 }
100
101 //____________________________________________________________________________
102 AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
103 {
104   // copy ctor
105   ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
106 }
107
108 //____________________________________________________________________________
109 void AliPHOSAnalyze::Copy(TObject & obj)
110 {
111   // copy an analysis into an other one
112   TObject::Copy(obj) ;
113   // I do nothing more because the copy is silly but the Code checkers requires one
114 }
115
116 //____________________________________________________________________________
117 AliPHOSAnalyze::~AliPHOSAnalyze()
118 {
119   // dtor
120
121   if(fRootFile->IsOpen()) fRootFile->Close() ; 
122   if(fRootFile)   {delete fRootFile ; fRootFile=0 ;}
123   if(fPHOS)       {delete fPHOS     ; fPHOS    =0 ;}
124   if(fClu)        {delete fClu      ; fClu     =0 ;}
125   if(fPID)        {delete fPID      ; fPID     =0 ;}
126   if(fRec)        {delete fRec      ; fRec     =0 ;}
127   if(fTrs)        {delete fTrs      ; fTrs     =0 ;}
128
129 }
130 //____________________________________________________________________________
131 void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){
132   //Draws pimary particles and reconstructed 
133   //digits, RecPoints, RecPartices etc 
134   //for event Nevent in the module Nmod.
135
136   TH2F * digitOccupancy  = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.);
137   TH2F * emcOccupancy    = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.);
138   TH2F * ppsdUp          = new TH2F("ppsdUp","PPSD Up digits",     128,-71.,71.,128,-71.,71.) ;
139   TH2F * ppsdUpCl        = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ;
140   TH2F * ppsdLow         = new TH2F("ppsdLow","PPSD Low digits",     128,-71.,71.,128,-71.,71.) ;
141   TH2F * ppsdLowCl       = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ;
142   TH2F * nbar            = new TH2F("nbar","Primary nbar",    64,-71.,71.,64,-71.,71.);
143   TH2F * phot            = new TH2F("phot","Primary Photon",  64,-71.,71.,64,-71.,71.);
144   TH2F * charg           = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.);
145   TH2F * recPhot         = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.);
146   TH2F * recNbar         = new TH2F("recNbar","RecParticles with primary Nbar",  64,-71.,71.,64,-71.,71.);
147
148   //========== Create the Clusterizer
149   fClu = new AliPHOSClusterizerv1() ; 
150
151   fClu->SetEmcEnergyThreshold(0.05) ;
152   fClu->SetEmcClusteringThreshold(0.20) ;
153   fClu->SetPpsdEnergyThreshold    (0.0000002) ;
154   fClu->SetPpsdClusteringThreshold(0.0000001) ;
155   fClu->SetLocalMaxCut(0.03) ;
156   fClu->SetCalibrationParameters(0., 0.00000001) ;
157   
158   gAlice->GetEvent(Nevent);
159   
160   TClonesArray * primaryList  = gAlice->Particles();
161   
162   TParticle * primary ;
163   Int_t iPrimary ;
164   for ( iPrimary = 0 ; iPrimary < primaryList->GetEntries() ; iPrimary++)
165     {
166       primary = (TParticle*)primaryList->At(iPrimary) ;
167       Int_t primaryType = primary->GetPdgCode() ;
168       if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
169         Int_t moduleNumber ;
170         Double_t primX, primZ ;
171         fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
172         if(moduleNumber==Nmod)
173           charg->Fill(primZ,primX,primary->Energy()) ;
174       }
175       if( primaryType == 22 ) {
176         Int_t moduleNumber ;
177         Double_t primX, primZ ;
178         fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
179         if(moduleNumber==Nmod)
180           phot->Fill(primZ,primX,primary->Energy()) ;
181       }
182       else{
183         if( primaryType == -2112 ) {
184           Int_t moduleNumber ;
185           Double_t primX, primZ ;
186           fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
187           if(moduleNumber==Nmod)
188             nbar->Fill(primZ,primX,primary->Energy()) ;
189         }
190       }
191     }  
192
193   fPHOS->SetTreeAddress() ;
194
195   gAlice->TreeD()->GetEvent(0) ;
196   gAlice->TreeR()->GetEvent(0) ;
197   
198   TObjArray ** emcRecPoints =  fPHOS->EmcRecPoints() ;
199   TObjArray ** ppsdRecPoints = fPHOS->PpsdRecPoints() ;
200   TClonesArray ** recParticleList  = fPHOS->RecParticles() ;
201   
202   Int_t iDigit ;
203   AliPHOSDigit * digit ;
204   
205   for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
206     {
207       digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
208       Int_t relid[4];
209       fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
210       Float_t x,z ;
211       fGeom->RelPosInModule(relid,x,z) ;
212       Float_t e = fClu->Calibrate(digit->GetAmp()) ;
213       if(relid[0]==Nmod){
214         if(relid[1]==0)  //EMC
215           digitOccupancy->Fill(x,z,e) ;
216         if((relid[1]>0)&&(relid[1]<17))
217           ppsdUp->Fill(x,z,e) ;
218         if(relid[1]>16)
219           ppsdLow->Fill(x,z,e) ;
220       }
221     }
222   
223   Int_t irecp ;
224   TVector3 pos ;
225   
226   for(irecp = 0; irecp < (*emcRecPoints)->GetEntries() ; irecp ++){
227     AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint*)(*emcRecPoints)->At(irecp) ;
228     if(emc->GetPHOSMod()==Nmod){
229       emc->GetLocalPosition(pos) ;
230       emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
231     }
232   }
233   
234   for(irecp = 0; irecp < (*ppsdRecPoints)->GetEntries() ; irecp ++){
235     AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)(*ppsdRecPoints)->At(irecp) ;
236     if(ppsd->GetPHOSMod()==Nmod){
237       ppsd->GetLocalPosition(pos) ;
238       if(ppsd->GetUp())
239         ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
240       else
241         ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
242     }
243   }
244   
245   AliPHOSRecParticle * recParticle ;
246   Int_t iRecParticle ;
247   for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
248     {
249       recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
250       
251       Int_t moduleNumberRec ;
252       Double_t recX, recZ ;
253       fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
254       if(moduleNumberRec == Nmod){
255         
256         Double_t minDistance = 5. ;
257         Int_t closestPrimary = -1 ;
258         
259         Int_t numberofprimaries ;
260         Int_t * listofprimaries  = recParticle->GetPrimaries(numberofprimaries)  ;
261         Int_t index ;
262         TParticle * primary ;
263         Double_t distance = minDistance ;
264           
265         for ( index = 0 ; index < numberofprimaries ; index++){
266           primary = (TParticle*)primaryList->At(listofprimaries[index]) ;
267           Int_t moduleNumber ;
268           Double_t primX, primZ ;
269           fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
270           if(moduleNumberRec == moduleNumber)
271             distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
272           if(minDistance > distance)
273             {
274               minDistance = distance ;
275               closestPrimary = listofprimaries[index] ;
276             }
277         }
278         
279         if(closestPrimary >=0 ){
280           
281           Int_t primaryType = ((TParticle *)primaryList->At(closestPrimary))->GetPdgCode() ;
282           
283           if(primaryType==22)
284             recPhot->Fill(recZ,recX,recParticle->Energy()) ;
285           else
286             if(primaryType==-2112)
287               recNbar->Fill(recZ,recX,recParticle->Energy()) ; 
288         }
289       }
290     }
291   
292   
293   digitOccupancy->Draw("box") ;
294   emcOccupancy->SetLineColor(2) ;
295   emcOccupancy->Draw("boxsame") ;
296   ppsdUp->SetLineColor(3) ;
297   ppsdUp->Draw("boxsame") ;
298   ppsdLow->SetLineColor(4) ;
299   ppsdLow->Draw("boxsame") ;
300   phot->SetLineColor(8) ;
301   phot->Draw("boxsame") ;
302   nbar->SetLineColor(6) ;
303   nbar->Draw("boxsame") ;
304   
305 }
306 //____________________________________________________________________________
307  void AliPHOSAnalyze::Reconstruct(Int_t nevents,Int_t firstEvent )    
308 {     
309
310   // Performs reconstruction of EMC and CPV (GPS2, IHEP or MIXT)
311   // for events from FirstEvent to Nevents
312
313   Int_t ievent ;   
314   for ( ievent=firstEvent; ievent<nevents; ievent++) {  
315     if (ievent==firstEvent) {
316       cout << "Analyze > Starting Reconstructing " << endl ; 
317       //========== Create the Clusterizer
318       fClu = new AliPHOSClusterizerv1() ; 
319       fClu->SetEmcEnergyThreshold(0.05) ; 
320       fClu->SetEmcClusteringThreshold(0.20) ; 
321       fClu->SetLocalMaxCut(0.03) ;
322       if      (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
323         fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
324         fClu->SetPpsdClusteringThreshold(0.0000001) ; 
325       }
326       else if (strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
327         fClu->SetLocalMaxCutCPV(0.03) ;
328         fClu->SetLogWeightCutCPV(4.0) ;
329         fClu->SetCpvEnergyThreshold(0.09) ;
330       }
331       fClu->SetCalibrationParameters(0., 0.00000001) ; 
332       
333       //========== Creates the track segment maker
334       fTrs = new AliPHOSTrackSegmentMakerv1()  ;
335           //      fTrs->UnsetUnfoldFlag() ; 
336      
337       //========== Creates the particle identifier for GPS2 only
338       if      (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
339         fPID = new AliPHOSPIDv1() ;
340         fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ; 
341       }   
342       
343       //========== Creates the Reconstructioner
344       fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; 
345       if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);     
346     }
347       
348     if (fDebugLevel != 0 ||
349         (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
350       cout <<  "======= Analyze ======> Event " << ievent+1 << endl ;
351     
352     //=========== Connects the various Tree's for evt
353     Int_t tracks = gAlice->GetEvent(ievent);
354
355     fPHOS->Hit2Digit(tracks) ;
356         
357     //=========== Do the reconstruction
358     fPHOS->Reconstruction(fRec);    
359
360   }
361   
362   if(fClu)      {delete fClu      ; fClu     =0 ;}
363   if(fPID)      {delete fPID      ; fPID     =0 ;}
364   if(fRec)      {delete fRec      ; fRec     =0 ;}
365   if(fTrs)      {delete fTrs      ; fTrs     =0 ;}
366   
367 }
368
369 //-------------------------------------------------------------------------------------
370 void AliPHOSAnalyze::ReadAndPrintCPV(Int_t EvFirst, Int_t EvLast)
371 {
372   //
373   // Read and print generated and reconstructed hits in CPV
374   // for events from EvFirst to Nevent.
375   // If only EvFirst is defined, print only this one event.
376   // Author: Yuri Kharlov
377   // 12 October 2000
378   //
379
380   if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
381   for ( Int_t ievent=EvFirst; ievent<=EvLast; ievent++) {  
382     
383     //========== Event Number>
384     cout << endl <<  "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
385     
386     //=========== Connects the various Tree's for evt
387     Int_t ntracks = gAlice->GetEvent(ievent);
388
389     //========== Creating branches ===================================
390     AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
391     gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints  ) ;
392     
393     AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->CpvRecPoints() ;
394     gAlice->TreeR()->SetBranchAddress( "PHOSCpvRP", cpvRecPoints ) ;
395
396     // Read and print CPV hits
397       
398     AliPHOSCPVModule cpvModule;
399     TClonesArray    *cpvHits;
400     Int_t           nCPVhits;
401     AliPHOSCPVHit   *cpvHit;
402     TLorentzVector   p;
403     Float_t          xgen, zgen;
404     Int_t            ipart;
405     Int_t            nGenHits = 0;
406     for (Int_t itrack=0; itrack<ntracks; itrack++) {
407       //=========== Get the Hits Tree for the Primary track itrack
408       gAlice->ResetHits();
409       gAlice->TreeH()->GetEvent(itrack);
410       Int_t iModule = 0 ;       
411       for (iModule=0; iModule < fGeom->GetNCPVModules(); iModule++) {
412         cpvModule = fPHOS->GetCPVModule(iModule);
413         cpvHits   = cpvModule.Hits();
414         nCPVhits  = cpvHits->GetEntriesFast();
415         for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
416           nGenHits++;
417           cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
418           p      = cpvHit->GetMomentum();
419           xgen   = cpvHit->X();
420           zgen   = cpvHit->Y();
421           ipart  = cpvHit->GetIpart();
422           printf("CPV hit in module %d: ",iModule+1);
423           printf(" p = (%f, %f, %f, %f) GeV,\n",
424                  p.Px(),p.Py(),p.Pz(),p.Energy());
425           printf("                  (X,Z) = (%8.4f, %8.4f) cm, ipart = %d\n",
426                  xgen,zgen,ipart);
427         }
428       }
429     }
430
431     // Read and print CPV reconstructed points
432
433     //=========== Gets the Reconstruction TTree
434     gAlice->TreeR()->GetEvent(0) ;
435     printf("Recpoints: %d\n",(*fPHOS->CpvRecPoints())->GetEntries());
436     TIter nextRP(*fPHOS->CpvRecPoints() ) ;
437     AliPHOSCpvRecPoint *cpvRecPoint ;
438     Int_t nRecPoints = 0;
439     while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
440       nRecPoints++;
441       TVector3  locpos;
442       cpvRecPoint->GetLocalPosition(locpos);
443       Int_t phosModule = cpvRecPoint->GetPHOSMod();
444       printf("CPV recpoint in module %d: (X,Z) = (%f,%f) cm\n",
445              phosModule,locpos.X(),locpos.Z());
446     }
447     printf("This event has %d generated hits and %d reconstructed points\n",
448            nGenHits,nRecPoints);
449   }
450 }
451
452 //____________________________________________________________________________
453 void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
454 {
455   //
456   // Analyzes CPV characteristics
457   // Author: Yuri Kharlov
458   // 9 October 2000
459   //
460
461   // Book histograms
462
463   TH1F *hDx   = new TH1F("hDx"  ,"CPV x-resolution@reconstruction",100,-5. , 5.);
464   TH1F *hDz   = new TH1F("hDz"  ,"CPV z-resolution@reconstruction",100,-5. , 5.);
465   TH1F *hDr   = new TH1F("hDr"  ,"CPV r-resolution@reconstruction",100, 0. , 5.);
466   TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",      21,-0.5,20.5);
467   TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length"  ,      21,-0.5,20.5);
468   TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length"    ,      21,-0.5,20.5);
469
470   cout << "Start CPV Analysis"<< endl ;
471   for ( Int_t ievent=0; ievent<Nevents; ievent++) {  
472       
473     //========== Event Number>         
474 //      if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
475       cout << endl <<  "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
476     
477     //=========== Connects the various Tree's for evt
478     Int_t ntracks = gAlice->GetEvent(ievent);
479     
480     //========== Creating branches ===================================
481     AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
482     gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints  ) ;
483     
484     AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
485     gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
486
487     // Create and fill arrays of hits for each CPV module
488       
489     Int_t nOfModules = fGeom->GetNModules();
490     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
491     Int_t iModule = 0;  
492     for (iModule=0; iModule < nOfModules; iModule++)
493       hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
494
495     AliPHOSCPVModule cpvModule;
496     TClonesArray    *cpvHits;
497     Int_t           nCPVhits;
498     AliPHOSCPVHit   *cpvHit;
499     TLorentzVector   p;
500     Float_t          xzgen[2];
501     Int_t            ipart;
502
503     // First go through all primary tracks and fill the arrays
504     // of hits per each CPV module
505
506     for (Int_t itrack=0; itrack<ntracks; itrack++) {
507       // Get the Hits Tree for the Primary track itrack
508       gAlice->ResetHits();
509       gAlice->TreeH()->GetEvent(itrack);
510       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
511         cpvModule = fPHOS->GetCPVModule(iModule);
512         cpvHits   = cpvModule.Hits();
513         nCPVhits  = cpvHits->GetEntriesFast();
514         for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
515           cpvHit   = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
516           p        = cpvHit->GetMomentum();
517           xzgen[0] = cpvHit->X();
518           xzgen[1] = cpvHit->Y();
519           ipart    = cpvHit->GetIpart();
520           TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
521           new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit);
522         }
523         cpvModule.Clear();
524       }
525     }
526     for (iModule=0; iModule < nOfModules; iModule++) {
527       Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
528       printf("Module %d has %d hits\n",iModule,nsum);
529     }
530
531     // Then go through reconstructed points and for each find
532     // the closeset hit
533     // The distance from the rec.point to the closest hit
534     // gives the coordinate resolution of the CPV
535
536     // Get the Reconstruction Tree
537     gAlice->TreeR()->GetEvent(0) ;
538     TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
539     AliPHOSCpvRecPoint *cpvRecPoint ;
540     Float_t xgen, zgen;
541     while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
542       TVector3  locpos;
543       cpvRecPoint->GetLocalPosition(locpos);
544       Int_t phosModule = cpvRecPoint->GetPHOSMod();
545       Int_t rpMult     = cpvRecPoint->GetDigitsMultiplicity();
546       Int_t rpMultX, rpMultZ;
547       cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
548       Float_t xrec  = locpos.X();
549       Float_t zrec  = locpos.Z();
550       Float_t dxmin = 1.e+10;
551       Float_t dzmin = 1.e+10;
552       Float_t r2min = 1.e+10;
553       Float_t r2;
554
555       cpvHits = hitsPerModule[phosModule-1];
556       Int_t nCPVhits  = cpvHits->GetEntriesFast();
557       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
558         cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
559         xgen   = cpvHit->X();
560         zgen   = cpvHit->Y();
561         r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2);
562         if ( r2 < r2min ) {
563           r2min = r2;
564           dxmin = xgen - xrec;
565           dzmin = zgen - zrec;
566         }
567       }
568       hDx  ->Fill(dxmin);
569       hDz  ->Fill(dzmin);
570       hDr  ->Fill(TMath::Sqrt(r2min));
571       hNrp ->Fill(rpMult);
572       hNrpX->Fill(rpMultX);
573       hNrpZ->Fill(rpMultZ);
574     }
575     delete [] hitsPerModule;
576   }
577   // Save histograms
578
579   Text_t outputname[80] ;
580   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
581   TFile output(outputname,"RECREATE");
582   output.cd();
583
584   hDx  ->Write() ;
585   hDz  ->Write() ;
586   hDr  ->Write() ;
587   hNrp ->Write() ;
588   hNrpX->Write() ;
589   hNrpZ->Write() ;
590
591   // Plot histograms
592
593   TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400);
594   gStyle->SetOptStat(111111);
595   gStyle->SetOptFit(1);
596   gStyle->SetOptDate(1);
597   cpvCanvas->Divide(3,2);
598
599   cpvCanvas->cd(1);
600   gPad->SetFillColor(10);
601   hNrp->SetFillColor(16);
602   hNrp->Draw();
603
604   cpvCanvas->cd(2);
605   gPad->SetFillColor(10);
606   hNrpX->SetFillColor(16);
607   hNrpX->Draw();
608
609   cpvCanvas->cd(3);
610   gPad->SetFillColor(10);
611   hNrpZ->SetFillColor(16);
612   hNrpZ->Draw();
613
614   cpvCanvas->cd(4);
615   gPad->SetFillColor(10);
616   hDx->SetFillColor(16);
617   hDx->Fit("gaus");
618   hDx->Draw();
619
620   cpvCanvas->cd(5);
621   gPad->SetFillColor(10);
622   hDz->SetFillColor(16);
623   hDz->Fit("gaus");
624   hDz->Draw();
625
626   cpvCanvas->cd(6);
627   gPad->SetFillColor(10);
628   hDr->SetFillColor(16);
629   hDr->Draw();
630
631   cpvCanvas->Print("CPV.ps");
632
633 }
634
635 //____________________________________________________________________________
636  void AliPHOSAnalyze::InvariantMass(Int_t Nevents )    
637 {
638   // Calculates Real and Mixed invariant mass distributions
639
640   const Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution 
641   Int_t mixedLoops = (Int_t )TMath::Ceil(Nevents/nMixedEvents) ;
642   
643   //========== Booking Histograms
644   TH2D * hRealEM   = new TH2D("hRealEM",   "Real for EM particles",      250,0.,1.,40,0.,4.) ;
645   TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
646   TH2D * hMixedEM  = new TH2D("hMixedEM",  "Mixed for EM particles",     250,0.,1.,40,0.,4.) ;
647   TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
648   
649   Int_t ievent;
650   Int_t eventInMixedLoop ;
651   
652   Int_t nRecParticles[4];//nMixedEvents] ;
653   
654   AliPHOSRecParticle::RecParticlesList * allRecParticleList  = new TClonesArray("AliPHOSRecParticle", nMixedEvents*1000) ;
655   
656   for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++  ){
657     Int_t iRecPhot = 0 ;
658     
659     for ( ievent=0; ievent < nMixedEvents; ievent++){        
660       
661       Int_t absEventNumber = eventInMixedLoop*nMixedEvents + ievent ;
662       
663       //=========== Connects the various Tree's for evt
664       gAlice->GetEvent(absEventNumber);
665
666       //========== Creating branches ===================================       
667       fPHOS->SetTreeAddress() ;
668       
669       gAlice->TreeD()->GetEvent(0) ;
670       gAlice->TreeR()->GetEvent(0) ;
671       
672       TClonesArray ** recParticleList  = fPHOS->RecParticles() ;
673       
674             
675       AliPHOSRecParticle * recParticle ;
676       Int_t iRecParticle ;
677       for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
678         {
679           recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
680           if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
681              (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){ 
682             new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ;
683             iRecPhot++;
684           }
685         }
686       
687         nRecParticles[ievent] = iRecPhot-1 ;  
688     }
689     
690     //Now calculate invariant mass:
691     Int_t irp1,irp2 ;
692     Int_t nCurEvent = 0 ;
693
694     for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
695       AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
696
697       for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
698         AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
699             
700         Double_t invMass ;
701         invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
702           (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
703           (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
704           (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
705         
706         if(invMass> 0)
707           invMass = TMath::Sqrt(invMass);
708         
709         Double_t pt ; 
710         pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
711
712         if(irp1 > nRecParticles[nCurEvent])
713           nCurEvent++;
714             
715         if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
716           hRealEM->Fill(invMass,pt);
717           if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
718             hRealPhot->Fill(invMass,pt);
719         }
720         else{
721           hMixedEM->Fill(invMass,pt);
722           if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
723             hMixedPhot->Fill(invMass,pt);
724         } //real-mixed
725             
726       } //loop over second rp
727     }//loop over first rp
728     allRecParticleList->Delete() ;
729   } //Loop over events
730   
731   delete allRecParticleList ;
732   
733   //writing output
734   TFile output("invmass.root","RECREATE");
735   output.cd();
736   
737   hRealEM->Write() ;
738   hRealPhot->Write() ;
739   hMixedEM->Write() ;
740   hMixedPhot->Write() ;
741   
742   output.Write();
743   output.Close();
744
745 }
746
747 //____________________________________________________________________________
748  void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst, Int_t EvLast)    
749 {
750   //
751   // Read and print generated and reconstructed hits in EMC
752   // for events from EvFirst to Nevent.
753   // If only EvFirst is defined, print only this one event.
754   // Author: Yuri Kharlov
755   // 24 November 2000
756   //
757
758   if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
759   Int_t ievent;
760   for (ievent=EvFirst; ievent<=EvLast; ievent++) {  
761     
762     //========== Event Number>
763     cout << endl <<  "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ;
764
765     //=========== Connects the various Tree's for evt
766     Int_t ntracks = gAlice->GetEvent(ievent);
767     fPHOS->SetTreeAddress() ;
768     
769     gAlice->TreeD()->GetEvent(0) ;
770     gAlice->TreeR()->GetEvent(0) ;
771
772     // Loop over reconstructed particles
773       
774     TClonesArray ** recParticleList  = fPHOS->RecParticles() ;     
775     AliPHOSRecParticle * recParticle ;
776     Int_t iRecParticle ;
777     Int_t *primList;
778     Int_t nPrimary;
779     for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) {
780       recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
781       Float_t recE = recParticle->Energy();
782       primList     = recParticle->GetPrimaries(nPrimary);
783       Int_t moduleNumberRec ;
784       Double_t recX, recZ ;
785       fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
786       printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n",
787              moduleNumberRec,recX,recZ,recE,*primList);
788     }
789
790     // Read and print EMC hits from EMCn branches
791       
792     AliPHOSCPVModule emcModule;
793     TClonesArray    *emcHits;
794     Int_t           nEMChits;
795     AliPHOSCPVHit   *emcHit;
796     TLorentzVector   p;
797     Float_t          xgen, zgen;
798     Int_t            ipart, primary;
799     Int_t            nGenHits = 0;
800     for (Int_t itrack=0; itrack<ntracks; itrack++) {
801       //=========== Get the Hits Tree for the Primary track itrack
802       gAlice->ResetHits();
803       gAlice->TreeH()->GetEvent(itrack);
804       Int_t iModule = 0 ;
805       for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
806         emcModule = fPHOS->GetEMCModule(iModule);
807         emcHits   = emcModule.Hits();
808         nEMChits  = emcHits->GetEntriesFast();
809         for (Int_t ihit=0; ihit<nEMChits; ihit++) {
810           nGenHits++;
811           emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
812           p      = emcHit->GetMomentum();
813           xgen   = emcHit->X();
814           zgen   = emcHit->Y();
815           ipart  = emcHit->GetIpart();
816           primary= emcHit->GetTrack();
817           printf("EMC hit A: module %d, ",iModule+1);
818           printf("    p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
819                  p.Px(),p.Py(),p.Pz(),p.Energy());
820           printf("                     (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
821                  xgen,zgen,ipart,primary);
822         }
823       }
824     }
825
826 //      // Read and print EMC hits from PHOS branch
827
828 //      for (Int_t itrack=0; itrack<ntracks; itrack++) {
829 //        //=========== Get the Hits Tree for the Primary track itrack
830 //        gAlice->ResetHits();
831 //        gAlice->TreeH()->GetEvent(itrack);
832 //        TClonesArray *hits = fPHOS->Hits();
833 //        AliPHOSHit   *hit ;
834 //        Int_t ihit;
835 //        for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
836 //      hit = (AliPHOSHit*)hits->At(ihit) ;
837 //      Float_t hitXYZ[3];
838 //      hitXYZ[0]   = hit->X();
839 //      hitXYZ[1]   = hit->Y();
840 //      hitXYZ[2]   = hit->Z();
841 //      ipart       = hit->GetPid();
842 //      primary     = hit->GetPrimary();
843 //      Int_t absId = hit->GetId();
844 //      Int_t relId[4];
845 //      fGeom->AbsToRelNumbering(absId, relId) ;
846 //      Int_t module = relId[0];
847 //      if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
848 //        printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
849 //               module,hitXYZ[0],hitXYZ[2],ipart,primary);
850 //        }
851 //      }
852
853   }
854 }
855
856 //____________________________________________________________________________
857  void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents)
858 {
859   //
860   // Read generated and reconstructed hits in EMC for Nevents events.
861   // Plots the coordinate and energy resolution histograms.
862   // Coordinate resolution is a difference between the reconstructed
863   // coordinate and the exact coordinate on the face of the PHOS
864   // Author: Yuri Kharlov
865   // 27 November 2000
866   //
867
868   // Book histograms
869
870   TH1F *hDx1   = new TH1F("hDx1"  ,"EMC x-resolution", 100,-5. , 5.);
871   TH1F *hDz1   = new TH1F("hDz1"  ,"EMC z-resolution", 100,-5. , 5.);
872   TH1F *hDE1   = new TH1F("hDE1"  ,"EMC E-resolution", 100,-2. , 2.);
873
874   TH2F *hDx2   = new TH2F("hDx2"  ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
875   TH2F *hDz2   = new TH2F("hDz2"  ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
876   TH2F *hDE2   = new TH2F("hDE2"  ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
877
878   cout << "Start EMC Analysis"<< endl ;
879   for (Int_t ievent=0; ievent<Nevents; ievent++) {  
880       
881     //========== Event Number>         
882     if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
883       cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
884     
885     //=========== Connects the various Tree's for evt
886     Int_t ntracks = gAlice->GetEvent(ievent);
887
888     fPHOS->SetTreeAddress() ;
889     
890     gAlice->TreeD()->GetEvent(0) ;
891     gAlice->TreeR()->GetEvent(0) ;
892
893     // Create and fill arrays of hits for each EMC module
894       
895     Int_t nOfModules = fGeom->GetNModules();
896     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
897     Int_t iModule;
898     for (iModule=0; iModule < nOfModules; iModule++)
899       hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
900
901     AliPHOSCPVModule emcModule;
902     TClonesArray    *emcHits;
903     Int_t           nEMChits;
904     AliPHOSCPVHit   *emcHit;
905
906     // First go through all primary tracks and fill the arrays
907     // of hits per each EMC module
908
909     for (Int_t itrack=0; itrack<ntracks; itrack++) {
910       // Get the Hits Tree for the Primary track itrack
911       gAlice->ResetHits();
912       gAlice->TreeH()->GetEvent(itrack);
913       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
914         emcModule = fPHOS->GetEMCModule(iModule);
915         emcHits   = emcModule.Hits();
916         nEMChits  = emcHits->GetEntriesFast();
917         for (Int_t ihit=0; ihit<nEMChits; ihit++) {
918           emcHit   = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
919           TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
920           new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
921         }
922         emcModule.Clear();
923       }
924     }
925
926     // Loop over reconstructed particles
927       
928     TClonesArray ** recParticleList  = fPHOS->RecParticles() ;     
929     AliPHOSRecParticle * recParticle ;
930     Int_t nEMCrecs = (*recParticleList)->GetEntries();
931     if (nEMCrecs == 1) {
932       recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
933       Float_t recE = recParticle->Energy();
934       Int_t phosModule;
935       Double_t recX, recZ ;
936       fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
937
938       // for this rec.point take the hit list in the same PHOS module
939
940       emcHits = hitsPerModule[phosModule-1];
941       Int_t nEMChits  = emcHits->GetEntriesFast();
942       if (nEMChits == 1) {
943         Float_t genX, genZ, genE;
944         for (Int_t ihit=0; ihit<nEMChits; ihit++) {
945           emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
946           genX   = emcHit->X();
947           genZ   = emcHit->Y();
948           genE   = emcHit->GetMomentum().E();
949         }
950         Float_t dx = recX - genX;
951         Float_t dz = recZ - genZ;
952         Float_t de = recE - genE;
953         hDx1  ->Fill(dx);
954         hDz1  ->Fill(dz);
955         hDE1  ->Fill(de);
956         hDx2  ->Fill(genE,dx);
957         hDz2  ->Fill(genE,dz);
958         hDE2  ->Fill(genE,recE);
959       }
960     }
961     delete [] hitsPerModule;
962   }
963   // Save histograms
964
965   Text_t outputname[80] ;
966   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
967   TFile output(outputname,"RECREATE");
968   output.cd();
969
970   hDx1  ->Write() ;
971   hDz1  ->Write() ;
972   hDE1  ->Write() ;
973   hDx2  ->Write() ;
974   hDz2  ->Write() ;
975   hDE2  ->Write() ;
976
977   // Plot histograms
978
979   TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
980   gStyle->SetOptStat(111111);
981   gStyle->SetOptFit(1);
982   gStyle->SetOptDate(1);
983   emcCanvas->Divide(3,1);
984
985   emcCanvas->cd(1);
986   gPad->SetFillColor(10);
987   hDx1->SetFillColor(16);
988   hDx1->Draw();
989
990   emcCanvas->cd(2);
991   gPad->SetFillColor(10);
992   hDz1->SetFillColor(16);
993   hDz1->Draw();
994
995   emcCanvas->cd(3);
996   gPad->SetFillColor(10);
997   hDE1->SetFillColor(16);
998   hDE1->Draw();
999
1000   emcCanvas->Print("EMC.ps");
1001
1002 }
1003
1004 //____________________________________________________________________________
1005  void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )    
1006 {
1007   // analyzes Nevents events and calculate Energy and Position resolution as well as
1008   // probaility of correct indentifiing of the incident particle
1009
1010   //========== Booking Histograms
1011   cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ; 
1012   BookResolutionHistograms();
1013
1014   Int_t counter[9][5] ;     
1015   Int_t i1,i2,totalInd = 0 ;
1016   for(i1 = 0; i1<9; i1++)
1017     for(i2 = 0; i2<5; i2++)
1018       counter[i1][i2] = 0 ;
1019   
1020   Int_t totalPrimary = 0 ;
1021   Int_t totalRecPart = 0 ;
1022   Int_t totalRPwithPrim = 0 ;
1023   Int_t ievent;
1024
1025   cout << "Start Analysing"<< endl ;
1026   for ( ievent=0; ievent<Nevents; ievent++)
1027     {  
1028       
1029       //========== Event Number>         
1030       //      if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
1031         cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;  
1032       
1033       //=========== Connects the various Tree's for evt
1034       gAlice->GetEvent(ievent);
1035
1036       //=========== Gets the Kine TTree
1037       gAlice->TreeK()->GetEvent(0) ;
1038       
1039       //=========== Gets the list of Primari Particles
1040       TClonesArray * primaryList  = gAlice->Particles();     
1041
1042       TParticle * primary ;
1043       Int_t iPrimary ;
1044       for ( iPrimary = 0 ; iPrimary < primaryList->GetEntries() ; iPrimary++)
1045         {
1046           primary = (TParticle*)primaryList->UncheckedAt(iPrimary) ;
1047           Int_t primaryType = primary->GetPdgCode() ;
1048           if( primaryType == 22 ) {
1049             Int_t moduleNumber ;
1050             Double_t primX, primZ ;
1051             fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1052             if(moduleNumber){
1053               fhPrimary->Fill(primary->Energy()) ;
1054               if(primary->Energy() > 0.3)
1055                 totalPrimary++ ;
1056             }
1057           } 
1058         }
1059       
1060       fPHOS->SetTreeAddress() ;
1061       
1062       gAlice->TreeD()->GetEvent(0) ;
1063       gAlice->TreeR()->GetEvent(0) ;
1064       
1065       TClonesArray ** recParticleList  = fPHOS->RecParticles() ;     
1066       
1067       AliPHOSRecParticle * recParticle ;
1068       Int_t iRecParticle ;
1069       for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
1070         {
1071           recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
1072           fhAllRP->Fill(CorrectEnergy(recParticle->Energy())) ;
1073           
1074           Int_t moduleNumberRec ;
1075           Double_t recX, recZ ;
1076           fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
1077           
1078           Double_t minDistance  = 100. ;
1079           Int_t closestPrimary = -1 ;
1080           
1081           Int_t numberofprimaries ;
1082           Int_t * listofprimaries  = recParticle->GetPrimaries(numberofprimaries)  ;
1083           Int_t index ;
1084           TParticle * primary ;
1085           Double_t distance = minDistance ;
1086           Double_t dX, dZ; 
1087           Double_t dXmin = 0.; 
1088           Double_t dZmin = 0. ;
1089           for ( index = 0 ; index < numberofprimaries ; index++){
1090             primary = (TParticle*)primaryList->UncheckedAt(listofprimaries[index]) ;
1091             Int_t moduleNumber ;
1092             Double_t primX, primZ ;
1093             fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1094             if(moduleNumberRec == moduleNumber) {
1095               dX = recX - primX;
1096               dZ = recZ - primZ;
1097               distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
1098               if(minDistance > distance) {
1099                 minDistance = distance ;
1100                 dXmin = dX;
1101                 dZmin = dZ;
1102                 closestPrimary = listofprimaries[index] ;
1103               }
1104             }
1105           }
1106           totalRecPart++ ;
1107
1108           if(closestPrimary >=0 ){
1109             totalRPwithPrim++;
1110             
1111             Int_t primaryType = ((TParticle *)primaryList->At(closestPrimary))->GetPdgCode() ;
1112 //          TParticlePDG* pDGparticle = ((TParticle *)primaryList->At(closestPrimary))->GetPDG();
1113 //          Double_t charge =  PDGparticle->Charge() ;
1114 //          if(charge)
1115 //            cout <<"Primary " <<primaryType << " E " << ((TParticle *)primaryList->At(closestPrimary))->Energy() << endl ;
1116             Int_t primaryCode ;
1117             switch(primaryType)
1118               {
1119               case 22:
1120                 primaryCode = 0;  //Photon
1121                 fhAllEnergy   ->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy()) ;
1122                 fhAllPosition ->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), minDistance) ;
1123                 fhAllPositionX->Fill(dXmin);
1124                 fhAllPositionZ->Fill(dZmin);
1125                 break;
1126               case 11 :
1127                 primaryCode = 1;  //Electron
1128                 break;
1129               case -11 :
1130                 primaryCode = 1;  //positron
1131                 break;
1132               case 321 :
1133                 primaryCode = 4;  //K+
1134                 break;
1135               case -321 :
1136                 primaryCode = 4;  //K-
1137                 break;
1138               case 310 :
1139                 primaryCode = 4;  //K0s
1140                 break;
1141               case 130 :
1142                 primaryCode = 4;  //K0l
1143                 break;
1144               case 211 :
1145                 primaryCode = 2;  //K0l
1146                 break;
1147               case -211 :
1148                 primaryCode = 2;  //K0l
1149                 break;
1150               case 2212 :
1151                 primaryCode = 2;  //K0l
1152                 break;
1153               case -2212 :
1154                 primaryCode = 2;  //K0l
1155                 break;
1156               default:
1157                 primaryCode = 3; //ELSE
1158                 break;
1159               }
1160             
1161             switch(recParticle->GetType())
1162               {
1163               case AliPHOSFastRecParticle::kGAMMA:
1164                 if(primaryType == 22){
1165                   fhPhotEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ; 
1166                   fhEMEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ; 
1167                   fhPPSDEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ; 
1168
1169                   fhPhotPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
1170                   fhEMPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
1171                   fhPPSDPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
1172
1173                   fhPhotReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1174                   fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1175                   fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1176
1177                   fhPhotPhot->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1178                 }
1179                 if(primaryType == 2112){ //neutron
1180                   fhNReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1181                   fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1182                   fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1183                 }
1184                 
1185                 if(primaryType == -2112){ //neutron ~
1186                   fhNBarReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1187                   fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1188                   fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1189                   
1190                 }
1191                 if(primaryCode == 2){
1192                   fhChargedReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1193                   fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1194                   fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1195                 }
1196                 
1197                 fhAllReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1198                 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1199                 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1200                 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1201                 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1202                 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1203                 counter[0][primaryCode]++;
1204                 break;
1205               case  AliPHOSFastRecParticle::kELECTRON:
1206                 if(primaryType == 22){ 
1207                   fhPhotElec->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1208                   fhEMEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ; 
1209                   fhEMPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
1210                   fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1211                   fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1212                 }         
1213                 if(primaryType == 2112){ //neutron
1214                   fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1215                   fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1216                 }
1217                 
1218                 if(primaryType == -2112){ //neutron ~
1219                   fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1220                   fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1221                   
1222                 }
1223                 if(primaryCode == 2){
1224                   fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1225                   fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1226                 }
1227                 
1228                 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1229                 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1230                 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1231                 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1232                 counter[1][primaryCode]++;
1233                 break;
1234               case  AliPHOSFastRecParticle::kNEUTRALHA:
1235                 if(primaryType == 22) 
1236                   fhPhotNeuH->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1237
1238                 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;           
1239                 counter[2][primaryCode]++;
1240                 break ;
1241               case  AliPHOSFastRecParticle::kNEUTRALEM:
1242                 if(primaryType == 22){
1243                   fhEMEnergy->Fill(((TParticle *)primaryList->At(closestPrimary))->Energy(),recParticle->Energy() ) ; 
1244                   fhEMPosition->Fill(((TParticle *)primaryList->At(closestPrimary))->Energy(),minDistance ) ;
1245                 
1246                   fhPhotNuEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1247                   fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1248                 }
1249                 if(primaryType == 2112) //neutron
1250                   fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1251                 
1252                 if(primaryType == -2112) //neutron ~
1253                   fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1254                 
1255                 if(primaryCode == 2)
1256                   fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1257                 
1258                 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1259                 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1260                 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1261
1262                 counter[3][primaryCode]++;
1263                 break ;
1264               case  AliPHOSFastRecParticle::kCHARGEDHA:
1265                 if(primaryType == 22) //photon
1266                   fhPhotChHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1267                 
1268                 counter[4][primaryCode]++ ;
1269                 break ;
1270               case  AliPHOSFastRecParticle::kGAMMAHA:
1271                   if(primaryType == 22){ //photon
1272                     fhPhotGaHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1273                     fhPPSDEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ; 
1274                     fhPPSDPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
1275                     fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1276                   }
1277                   if(primaryType == 2112){ //neutron
1278                     fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1279                   }
1280                 
1281                   if(primaryType == -2112){ //neutron ~
1282                     fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; 
1283                   }
1284                   if(primaryCode == 2){
1285                     fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1286                   }
1287                 
1288                   fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1289                   fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1290                   fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1291                   counter[5][primaryCode]++ ;
1292                   break ;       
1293               case  AliPHOSFastRecParticle::kABSURDEM:        
1294                 counter[6][primaryCode]++ ;
1295                 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1296                 break;
1297               case  AliPHOSFastRecParticle::kABSURDHA:
1298                 counter[7][primaryCode]++ ;
1299                 break;
1300               default:
1301                 counter[8][primaryCode]++ ;
1302                 break;
1303               }
1304           }
1305         }  
1306     }   // endfor
1307   SaveHistograms();
1308   cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
1309   cout << "Resolutions: Total primary       " << totalPrimary << endl ;
1310   cout << "Resoluitons: Total reconstracted " << totalRecPart << endl ;
1311   cout << "TotalReconstructed with Primarie " << totalRPwithPrim << endl ;
1312   cout << "                        Primary:   Photon   Electron   Ch. Hadr.  Neutr. Hadr  Kaons" << endl ; 
1313   cout << "             Detected as photon       " << counter[0][0] << "          " << counter[0][1] << "          " << counter[0][2] << "          " <<counter[0][3] << "          " << counter[0][4] << endl ;
1314   cout << "           Detected as electron       " << counter[1][0] << "          " << counter[1][1] << "          " << counter[1][2] << "          " <<counter[1][3] << "          " << counter[1][4] << endl ; 
1315   cout << "     Detected as neutral hadron       " << counter[2][0] << "          " << counter[2][1] << "          " << counter[2][2] << "          " <<counter[2][3] << "          " << counter[2][4] << endl ;
1316   cout << "         Detected as neutral EM       " << counter[3][0] << "          " << counter[3][1] << "          " << counter[3][2] << "          " <<counter[3][3] << "          " << counter[3][4] << endl ;
1317   cout << "     Detected as charged hadron       " << counter[4][0] << "          " << counter[4][1] << "          " << counter[4][2] << "          " <<counter[4][3] << "          " << counter[4][4] << endl ;
1318   cout << "       Detected as gamma-hadron       " << counter[5][0] << "          " << counter[5][1] << "          " << counter[5][2] << "          " <<counter[5][3] << "          " << counter[5][4] << endl ;
1319   cout << "          Detected as Absurd EM       " << counter[6][0] << "          " << counter[6][1] << "          " << counter[6][2] << "          " <<counter[6][3] << "          " << counter[6][4] << endl ;
1320   cout << "      Detected as absurd hadron       " << counter[7][0] << "          " << counter[7][1] << "          " << counter[7][2] << "          " <<counter[7][3] << "          " << counter[7][4] << endl ;
1321   cout << "          Detected as undefined       " << counter[8][0] << "          " << counter[8][1] << "          " << counter[8][2] << "          " <<counter[8][3] << "          " << counter[8][4] << endl ;
1322       
1323       for(i1 = 0; i1<9; i1++)
1324         for(i2 = 0; i2<5; i2++)
1325           totalInd+=counter[i1][i2] ;
1326       cout << "Indentified particles            " << totalInd << endl ;
1327       
1328 }           // endfunction
1329
1330
1331 //____________________________________________________________________________
1332 void  AliPHOSAnalyze::BookingHistograms()
1333 {
1334   // Books the histograms where the results of the analysis are stored (to be changed)
1335
1336   delete fhEmcDigit  ;
1337   delete fhVetoDigit  ;
1338   delete fhConvertorDigit   ;
1339   delete  fhEmcCluster   ;
1340   delete fhVetoCluster   ;
1341   delete fhConvertorCluster  ;
1342   delete fhConvertorEmc  ;
1343   
1344   fhEmcDigit                = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
1345   fhVetoDigit               = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
1346   fhConvertorDigit          = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
1347   fhEmcCluster              = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
1348   fhVetoCluster             = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
1349   fhConvertorCluster        = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
1350   fhConvertorEmc            = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
1351
1352 }
1353 //____________________________________________________________________________
1354 void  AliPHOSAnalyze::BookResolutionHistograms()
1355 {
1356   // Books the histograms where the results of the Resolution analysis are stored
1357
1358 //   if(fhAllEnergy)
1359 //     delete fhAllEnergy ;
1360 //   if(fhPhotEnergy)
1361 //     delete fhPhotEnergy ;
1362 //   if(fhEMEnergy)
1363 //     delete fhEMEnergy ;
1364 //   if(fhPPSDEnergy)
1365 //     delete fhPPSDEnergy ;
1366
1367
1368   fhAllEnergy  = new TH2F("hAllEnergy",  "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1369   fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1370   fhEMEnergy   = new TH2F("hEMEnergy",   "Energy of EM with primary photon",    100, 0., 5., 100, 0., 5.);
1371   fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon",  100, 0., 5., 100, 0., 5.);
1372
1373 //   if(fhAllPosition)
1374 //     delete fhAllPosition ;
1375 //   if(fhPhotPosition)
1376 //     delete fhPhotPosition ;
1377 //   if(fhEMPosition)
1378 //     delete fhEMPosition ;
1379 //   if(fhPPSDPosition)
1380 //     delete fhPPSDPosition ;
1381
1382
1383   fhAllPosition  = new TH2F("hAllPosition",  "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1384   fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1385   fhEMPosition   = new TH2F("hEMPosition",   "Position of EM with primary photon",    100, 0., 5., 100, 0., 5.);
1386   fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon",  100, 0., 5., 100, 0., 5.);
1387
1388   fhAllPositionX = new TH1F("hAllPositionX", "#Delta X of any RP with primary photon",100, -2., 2.);
1389   fhAllPositionZ = new TH1F("hAllPositionZ", "#Delta X of any RP with primary photon",100, -2., 2.);
1390
1391 //   if(fhAllReg)
1392 //     delete fhAllReg ;
1393 //   if(fhPhotReg)
1394 //     delete fhPhotReg ;
1395 //   if(fhNReg)
1396 //     delete fhNReg ;
1397 //   if(fhNBarReg)
1398 //     delete fhNBarReg ;
1399 //   if(fhChargedReg)
1400 //     delete fhChargedReg ;
1401   
1402   fhAllReg    = new TH1F("hAllReg",    "All primaries registered as photon",  100, 0., 5.);
1403   fhPhotReg   = new TH1F("hPhotReg",   "Photon registered as photon",         100, 0., 5.);
1404   fhNReg      = new TH1F("hNReg",      "N registered as photon",              100, 0., 5.);
1405   fhNBarReg   = new TH1F("hNBarReg",   "NBar registered as photon",           100, 0., 5.);
1406   fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
1407   
1408 //   if(fhAllEM)
1409 //     delete fhAllEM ;
1410 //   if(fhPhotEM)
1411 //     delete fhPhotEM ;
1412 //   if(fhNEM)
1413 //     delete fhNEM ;
1414 //   if(fhNBarEM)
1415 //     delete fhNBarEM ;
1416 //   if(fhChargedEM)
1417 //     delete fhChargedEM ;
1418   
1419   fhAllEM    = new TH1F("hAllEM",    "All primary registered as EM",100, 0., 5.);
1420   fhPhotEM   = new TH1F("hPhotEM",   "Photon registered as EM", 100, 0., 5.);
1421   fhNEM      = new TH1F("hNEM",      "N registered as EM",      100, 0., 5.);
1422   fhNBarEM   = new TH1F("hNBarEM",   "NBar registered as EM",   100, 0., 5.);
1423   fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
1424
1425 //   if(fhAllPPSD)
1426 //     delete fhAllPPSD ;
1427 //   if(fhPhotPPSD)
1428 //     delete fhPhotPPSD ;
1429 //   if(fhNPPSD)
1430 //     delete fhNPPSD ;
1431 //   if(fhNBarPPSD)
1432 //     delete fhNBarPPSD ;
1433 //   if(fhChargedPPSD)
1434 //     delete fhChargedPPSD ;
1435   
1436   fhAllPPSD    = new TH1F("hAllPPSD",    "All primary registered as PPSD",100, 0., 5.);
1437   fhPhotPPSD   = new TH1F("hPhotPPSD",   "Photon registered as PPSD", 100, 0., 5.);
1438   fhNPPSD      = new TH1F("hNPPSD",      "N registered as PPSD",      100, 0., 5.);
1439   fhNBarPPSD   = new TH1F("hNBarPPSD",   "NBar registered as PPSD",   100, 0., 5.);
1440   fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
1441   
1442 //   if(fhPrimary)
1443 //     delete fhPrimary ;
1444   fhPrimary= new TH1F("hPrimary", "hPrimary",  100, 0., 5.);
1445
1446 //   if(fhAllRP)
1447 //     delete fhAllRP ;
1448 //   if(fhVeto)
1449 //     delete fhVeto ;
1450 //   if(fhShape)
1451 //     delete fhShape ;
1452 //   if(fhPPSD)
1453 //     delete fhPPSD ;
1454
1455   fhAllRP = new TH1F("hAllRP","All Reconstructed particles",  100, 0., 5.);
1456   fhVeto  = new TH1F("hVeto", "All uncharged particles",      100, 0., 5.);
1457   fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.);
1458   fhPPSD  = new TH1F("hPPSD", "All PPSD photon particles",    100, 0., 5.);
1459
1460
1461 //   if(fhPhotPhot)
1462 //     delete fhPhotPhot ;
1463 //   if(fhPhotElec)
1464 //     delete fhPhotElec ;
1465 //   if(fhPhotNeuH)
1466 //     delete fhPhotNeuH ;
1467 //   if(fhPhotNuEM)
1468 //     delete fhPhotNuEM ;
1469 //   if(fhPhotChHa)
1470 //     delete fhPhotChHa ;
1471 //   if(fhPhotGaHa)
1472 //     delete fhPhotGaHa ;
1473
1474   fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.);   //Photon registered as photon
1475   fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.);   //Photon registered as Electron
1476   fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.);   //Photon registered as Neutral Hadron
1477   fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.);   //Photon registered as Neutral EM
1478   fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.);   //Photon registered as Charged Hadron
1479   fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.);   //Photon registered as Gamma-Hadron
1480 }
1481
1482 //____________________________________________________________________________
1483 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1484 {
1485   // Open the root file named "name"
1486   
1487   fRootFile   = new TFile(name, "update") ;
1488   return  fRootFile->IsOpen() ; 
1489 }
1490
1491 //____________________________________________________________________________
1492 void AliPHOSAnalyze::SaveHistograms()
1493 {
1494   // Saves the histograms in a root file named "name.analyzed" 
1495
1496   Text_t outputname[80] ;
1497   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1498   TFile output(outputname,"RECREATE");
1499   output.cd();
1500
1501   if (fhAllEnergy)    
1502     fhAllEnergy->Write() ;
1503   if (fhPhotEnergy)    
1504     fhPhotEnergy->Write() ;
1505   if(fhEMEnergy)
1506     fhEMEnergy->Write()  ;
1507   if(fhPPSDEnergy)
1508     fhPPSDEnergy->Write() ;
1509   if(fhAllPosition)
1510     fhAllPosition->Write() ;
1511   if(fhAllPositionX)
1512     fhAllPositionX->Write() ;
1513   if(fhAllPositionZ)
1514     fhAllPositionZ->Write() ;
1515   if(fhPhotPosition)
1516     fhPhotPosition->Write() ;
1517   if(fhEMPosition)
1518     fhEMPosition->Write() ;
1519   if(fhPPSDPosition)
1520     fhPPSDPosition->Write() ;
1521   if (fhAllReg) 
1522     fhAllReg->Write() ;
1523   if (fhPhotReg) 
1524     fhPhotReg->Write() ;
1525   if(fhNReg)
1526     fhNReg->Write() ;
1527   if(fhNBarReg)
1528     fhNBarReg->Write() ;
1529   if(fhChargedReg)
1530     fhChargedReg->Write() ;
1531   if (fhAllEM) 
1532     fhAllEM->Write() ;
1533   if (fhPhotEM) 
1534     fhPhotEM->Write() ;
1535   if(fhNEM)
1536     fhNEM->Write() ;
1537   if(fhNBarEM)
1538     fhNBarEM->Write() ;
1539   if(fhChargedEM)
1540     fhChargedEM->Write() ;
1541   if (fhAllPPSD) 
1542     fhAllPPSD->Write() ;
1543   if (fhPhotPPSD) 
1544     fhPhotPPSD->Write() ;
1545   if(fhNPPSD)
1546     fhNPPSD->Write() ;
1547   if(fhNBarPPSD)
1548     fhNBarPPSD->Write() ;
1549   if(fhChargedPPSD)
1550     fhChargedPPSD->Write() ;
1551   if(fhPrimary)
1552     fhPrimary->Write() ;
1553   if(fhAllRP)
1554     fhAllRP->Write()  ;
1555   if(fhVeto)
1556     fhVeto->Write()  ;
1557   if(fhShape)
1558     fhShape->Write()  ;
1559   if(fhPPSD)
1560     fhPPSD->Write()  ;
1561   if(fhPhotPhot)
1562     fhPhotPhot->Write() ;
1563   if(fhPhotElec)
1564     fhPhotElec->Write() ;
1565   if(fhPhotNeuH)
1566     fhPhotNeuH->Write() ;
1567   if(fhPhotNuEM)
1568     fhPhotNuEM->Write() ;
1569   if(fhPhotNuEM)
1570     fhPhotNuEM->Write() ;
1571   if(fhPhotChHa)
1572     fhPhotChHa->Write() ;
1573   if(fhPhotGaHa)
1574     fhPhotGaHa->Write() ;
1575   if(fhEnergyCorrelations)
1576     fhEnergyCorrelations->Write() ;
1577   
1578   output.Write();
1579   output.Close();
1580 }
1581 //____________________________________________________________________________
1582 Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
1583 {
1584   return ERecPart/0.8783 ;
1585 }
1586
1587 //____________________________________________________________________________
1588 void AliPHOSAnalyze::ResetHistograms()
1589 {
1590    fhEnergyCorrelations = 0 ;     //Energy correlations between Eloss in Convertor and PPSD(2)
1591
1592    fhEmcDigit = 0 ;               // Histo of digit energies in the Emc 
1593    fhVetoDigit = 0 ;              // Histo of digit energies in the Veto 
1594    fhConvertorDigit = 0 ;         // Histo of digit energies in the Convertor
1595    fhEmcCluster = 0 ;             // Histo of Cluster energies in Emc
1596    fhVetoCluster = 0 ;            // Histo of Cluster energies in Veto
1597    fhConvertorCluster = 0 ;       // Histo of Cluster energies in Convertor
1598    fhConvertorEmc = 0 ;           // 2d Convertor versus Emc energies
1599
1600    fhAllEnergy = 0 ;       
1601    fhPhotEnergy = 0 ;        // Total spectrum of detected photons
1602    fhEMEnergy = 0 ;         // Spectrum of detected electrons with electron primary
1603    fhPPSDEnergy = 0 ;
1604    fhAllPosition = 0 ; 
1605    fhAllPositionX = 0 ; 
1606    fhAllPositionZ = 0 ; 
1607    fhPhotPosition = 0 ; 
1608    fhEMPosition = 0 ; 
1609    fhPPSDPosition = 0 ; 
1610
1611    fhPhotReg = 0 ;          
1612    fhAllReg = 0 ;          
1613    fhNReg = 0 ;          
1614    fhNBarReg = 0 ;          
1615    fhChargedReg = 0 ;          
1616    fhPhotEM = 0 ;          
1617    fhAllEM = 0 ;          
1618    fhNEM = 0 ;          
1619    fhNBarEM = 0 ;          
1620    fhChargedEM = 0 ;          
1621    fhPhotPPSD = 0 ;          
1622    fhAllPPSD = 0 ;          
1623    fhNPPSD = 0 ;          
1624    fhNBarPPSD = 0 ;          
1625    fhChargedPPSD = 0 ;          
1626
1627    fhPrimary = 0 ;          
1628
1629    fhPhotPhot = 0 ;
1630    fhPhotElec = 0 ;
1631    fhPhotNeuH = 0 ;
1632    fhPhotNuEM = 0 ; 
1633    fhPhotChHa = 0 ;
1634    fhPhotGaHa = 0 ;
1635
1636 }