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