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