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