Transition to NewIO
[u/mrichter/AliRoot.git] / ZDC / AliZDC.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 //                                                                           //
20 //                      Zero Degree Calorimeter                              //
21 //           This class contains the basic functions for the ZDCs;           //
22 //            functions specific to one particular geometry are              //
23 //                      contained in the derived classes                     //
24 //                                                                           //
25 ///////////////////////////////////////////////////////////////////////////////
26
27 #include <stdlib.h>
28 #include <Riostream.h>
29
30 // --- ROOT system
31 #include <TBRIK.h>
32 #include <TDirectory.h>
33 #include <TF1.h>
34 #include <TFile.h>
35 #include <TGeometry.h>
36 #include <TNode.h>
37 #include <TTree.h>
38 #include <TVirtualMC.h>
39
40 // --- AliRoot header files
41 #include "AliDetector.h"
42 #include "AliZDC.h"
43 #include "AliZDCDigit.h"
44 #include "AliZDCHit.h"
45 #include "AliZDCMergedHit.h"
46 #include "AliZDCMerger.h"
47 #include "AliZDCReco.h"
48
49 #include "AliConst.h"
50
51 #include "AliHeader.h"
52 #include "AliLoader.h"
53 #include "AliRun.h"
54
55  
56 ClassImp(AliZDC)
57  
58 //_____________________________________________________________________________
59 AliZDC::AliZDC()
60 {
61   //
62   // Default constructor for the Zero Degree Calorimeter base class
63   //
64   
65   fIshunt     = 1;
66   fNoShower   = 0;
67
68   fMerger     = 0;
69   fHits       = 0;
70   fNhits      = 0;
71
72   fDigits     = 0;
73   fNdigits    = 0;
74
75   fMergedHits = 0;
76
77   fNRecPoints = 0;
78   fRecPoints  = 0;
79   
80 }
81  
82 //_____________________________________________________________________________
83 AliZDC::AliZDC(const char *name, const char *title)
84   : AliDetector(name,title)
85 {
86   //
87   // Standard constructor for the Zero Degree Calorimeter base class
88   //
89
90   fIshunt   = 1;
91   fNoShower = 0;
92   fMerger   = 0;
93
94   // Allocate the hits array  
95   fHits   = new TClonesArray("AliZDCHit",1000);
96   gAlice->AddHitList(fHits);
97   // Allocate the merged hits array  
98   fMergedHits = new TClonesArray("AliZDCMergedHit",1000);
99
100   // Allocate the digits array  
101   fDigits = new TClonesArray("AliZDCDigit",1000);
102   
103   fNRecPoints = 0;
104   fRecPoints = 0;  
105
106 }
107 //____________________________________________________________________________ 
108 AliZDC::~AliZDC()
109 {
110   //
111   // ZDC destructor
112   //
113
114   fIshunt   = 0;
115   
116   if(fMerger) delete fMerger;
117
118 }
119 //_____________________________________________________________________________
120 void AliZDC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
121 {
122   //
123   //            Add a ZDC hit to the hit list.
124   // -> We make use of 2 array of hits:
125   // [1]  fHits (the usual one) that contains hits for each PRIMARY
126   // [2]  fStHits that contains hits for each EVENT and is used to
127   //      obtain digits at the end of each event
128   //
129   
130   static Float_t primKinEn, xImpact, yImpact, sFlag;
131
132   AliZDCHit *newquad, *curprimquad;
133   newquad = new AliZDCHit(fIshunt, track, vol, hits);
134   TClonesArray &lhits = *fHits;
135   
136   if(fNhits==0){
137       // First hit -> setting flag for primary or secondary particle
138       Int_t primary = gAlice->GetPrimary(track);     
139       if(track != primary){
140         newquad->fSFlag = 1;  // SECONDARY particle entering the ZDC
141       }
142       else if(track == primary){
143         newquad->fSFlag = 0;  // PRIMARY particle entering the ZDC
144       }  
145       sFlag     = newquad->fSFlag;
146       primKinEn = newquad->fPrimKinEn;
147       xImpact   = newquad->fXImpact;
148       yImpact   = newquad->fYImpact;
149    }
150    else{       
151       newquad->fPrimKinEn = primKinEn;
152       newquad->fXImpact = xImpact;
153       newquad->fYImpact = yImpact;
154       newquad->fSFlag   = sFlag;
155    }
156  
157   Int_t j;
158   for(j=0; j<fNhits; j++){
159     // If hits are equal (same track, same volume), sum them.
160      curprimquad = (AliZDCHit*) lhits[j];
161      if(*curprimquad == *newquad){
162         *curprimquad = *curprimquad+*newquad;
163         delete newquad;
164         return;
165      } 
166   }
167
168     //Otherwise create a new hit
169     new(lhits[fNhits]) AliZDCHit(newquad);
170     fNhits++;
171     
172     delete newquad;
173 }
174
175 //_____________________________________________________________________________
176 void  AliZDC::AddDigit(Int_t *sect, Int_t digit)
177 {
178 //
179   AliZDCDigit *newdigit;
180   newdigit = new AliZDCDigit(sect, digit);
181   
182 //  printf("\n  AddDigit -> sector[0] = %d, sector[1] = %d, digit = %d",
183 //         sect[0], sect[1], digit);
184   new((*fDigits)[fNdigits]) AliZDCDigit(*newdigit);
185   fNdigits++;
186   delete newdigit;
187 }
188       
189 //_____________________________________________________________________________
190 void AliZDC::BuildGeometry()
191 {
192   //
193   // Build the ROOT TNode geometry for event display 
194   // in the Zero Degree Calorimeter
195   // This routine is dummy for the moment
196   //
197
198   TNode *node, *top;
199   TBRIK *brik;
200   const int kColorZDC  = kBlue;
201   
202   //
203   top=gAlice->GetGeometry()->GetNode("alice");
204   
205   // ZDC
206     brik = new TBRIK("S_ZDC","ZDC box","void",300,300,5);
207     top->cd();
208     node = new TNode("ZDC","ZDC","S_ZDC",0,0,600,"");
209     node->SetLineColor(kColorZDC);
210     fNodes->Add(node);
211 }
212
213 //_____________________________________________________________________________
214 Int_t AliZDC::DistancetoPrimitive(Int_t , Int_t )
215 {
216   //
217   // Distance from the mouse to the Zero Degree Calorimeter
218   // Dummy routine
219   //
220   return 9999;
221 }
222
223 //____________________________________________________________________________
224 Float_t AliZDC::ZMin(void) const
225 {
226   // Minimum dimension of the ZDC module in z
227   return 11600.;
228 }
229
230 //____________________________________________________________________________
231 Float_t AliZDC::ZMax(void) const
232 {
233   // Maximum dimension of the ZDC module in z
234   return  11750.;
235 }
236   
237
238 //_____________________________________________________________________________
239  void AliZDC::MakeBranch(Option_t *opt, const char *file)
240 {
241   //
242   // Create Tree branches for the ZDC
243   //
244
245   char branchname[10];
246   sprintf(branchname,"%s",GetName());
247
248   const char *cH = strstr(opt,"H");
249   
250   if (cH && fLoader->TreeH())
251    fHits   = new TClonesArray("AliZDCHit",1000); 
252   
253   AliDetector::MakeBranch(opt);
254
255   const char *cS = strstr(opt,"S");
256
257   if (fLoader->TreeS() && cS) {
258     if(fMergedHits!=0) fMergedHits->Clear();
259     else fMergedHits = new TClonesArray ("AliZDCMergedHit",1000);
260     MakeBranchInTree(fLoader->TreeS(), 
261                      branchname, &fMergedHits, fBufferSize, file) ;
262     printf("* AliZDC::MakeBranch    * Making Branch %s for SDigits\n\n",branchname);
263   }
264
265     
266   const char *cD = strstr(opt,"D");
267
268   if (fLoader->TreeD() && cD) {
269     if(fDigits!=0) fDigits->Clear();
270     else fDigits = new TClonesArray ("AliZDCDigit",1000);
271     MakeBranchInTree(fLoader->TreeD(), 
272                      branchname, &fDigits, fBufferSize, file) ;
273     printf("* AliZDC::MakeBranch    * Making Branch %s for Digits\n\n",branchname);
274   }
275
276   
277   const char *cR = strstr(opt,"R");
278
279   if (gAlice->TreeR() && cR) {
280     if(fRecPoints==0) fRecPoints = new TClonesArray("AliZDCReco",1000);
281     MakeBranchInTree(gAlice->TreeR(), 
282                      branchname, &fRecPoints, fBufferSize, file) ;
283     printf("* AliZDC::MakeBranch    * Making Branch %s for RecPoints\n\n",branchname);   }
284           
285 }
286
287 //_____________________________________________________________________________
288  void AliZDC::MakeBranchInTreeS(TTree *treeS, const char *file)
289 {
290   // MakeBranchInTree
291   const Int_t kBufferSize = 4000;
292   char  branchname[20];
293   sprintf(branchname,"%s",GetName());
294   if (fMergedHits==0x0) fMergedHits = new TClonesArray("AliZDCMergedHit",1000);
295   MakeBranchInTree(treeS, branchname, &fMergedHits, kBufferSize, file) ;
296   printf("* AliZDC::MakeBranch    * Making Branch %s for SDigits\n\n",branchname);
297
298 }
299 //_____________________________________________________________________________
300  void AliZDC::MakeBranchInTreeD(TTree *treeD, const char *file)
301 {
302   // MakeBranchInTree
303   const Int_t kBufferSize = 4000;
304   char  branchname[20];
305   sprintf(branchname,"%s",GetName());
306   if (fDigits == 0x0) fDigits = new TClonesArray("AliZDCDigit",1000);
307   MakeBranchInTree(treeD, branchname, &fDigits, kBufferSize, file) ;
308   printf("* AliZDC::MakeBranch    * Making Branch %s for Digits\n\n",branchname);
309
310 }
311 //_____________________________________________________________________________
312  void AliZDC::MakeBranchInTreeR(TTree *treeR, const char *file)
313 {
314   // MakeBranchInTree
315   const Int_t kBufferSize = 4000;
316   char  branchname[20];
317   sprintf(branchname,"%s",GetName());
318   MakeBranchInTree(treeR, branchname, &fRecPoints, kBufferSize, file) ;
319   printf("* AliZDC::MakeBranch    * Making Branch %s for RecPoints\n\n",branchname);
320
321 }
322 //_____________________________________________________________________________
323 void AliZDC::Hits2SDigits()
324 {
325   printf("\n    Entering AliZDC::SDigits2Digits() ");
326   
327   //----------------------------------------------------------------
328   if(!fMerger){ 
329     printf("    ZDC digitization (without merging)\n");
330
331     AliZDCMergedHit *MHit;
332     Int_t j, sector[2];
333     Float_t MHits[7];
334     fNMergedhits = 0;
335
336     TTree *treeH = TreeH();
337     Int_t ntracks = (Int_t) treeH->GetEntries();
338     gAlice->ResetHits();
339   
340     // Tracks loop
341     for(Int_t itrack=0; itrack<ntracks; itrack++){
342        treeH->GetEvent(itrack);
343        for(AliZDCHit* zdcHit=(AliZDCHit*)this->FirstHit(-1); zdcHit;
344                       zdcHit = (AliZDCHit*)this->NextHit()){ 
345                       
346            for(j=0; j<2; j++) sector[j] = zdcHit->GetVolume(j);
347            MHits[0] = zdcHit->GetPrimKinEn();
348            MHits[1] = zdcHit->GetXImpact();
349            MHits[2] = zdcHit->GetYImpact();
350            MHits[3] = zdcHit->GetSFlag();
351            MHits[4] = zdcHit->GetLightPMQ();
352            MHits[5] = zdcHit->GetLightPMC();
353            MHits[6] = zdcHit->GetEnergy();
354        }//Hits loop
355        
356           MHit = new AliZDCMergedHit(sector, MHits);
357           new((*fMergedHits)[fNMergedhits]) AliZDCMergedHit(*MHit);       
358           TClonesArray &sdigits = *fMergedHits;
359           new (sdigits[fNMergedhits]) AliZDCMergedHit(*MHit);
360           fNMergedhits++;
361           delete MHit;
362     }
363     gAlice->TreeS()->Fill();
364     gAlice->TreeS()->AutoSave(); 
365     gAlice->TreeS()->Reset();  
366   }
367   //----------------------------------------------------------------
368   else if(fMerger){
369     printf("    ZDC merging and digitization\n");
370     // ### Initialise merging
371     fMerger -> InitMerging();
372
373     // SDigits tree
374
375
376
377     TTree *treeS = fLoader->TreeS();
378     if (treeS == 0x0)
379      {
380       Int_t retval = fLoader->LoadSDigits();
381       if (retval)
382        {
383          Error("Hits2SDigits","Error while loading S. Digits");
384          return;
385        }
386       treeS = fLoader->TreeS();
387      }
388     
389     if(!treeS){
390       printf("\n ERROR -> Can't find TreeS%d in background file\n",fMerger->EvNum());
391     }    
392
393     // ### Get TCA of MergedHits from AliZDCMerger
394     fMergedHits  = fMerger->MergedHits();
395     fNMergedhits = fMerger->GetNMhits();
396
397     // Branch address
398     char branchSDname[20];
399     sprintf(branchSDname,"%s",GetName());
400     if(treeS && fMergedHits){
401       TBranch *branchSD = treeS->GetBranch(branchSDname);
402       if(branchSD) branchSD->SetAddress(&fMergedHits);
403       else if(!branchSD) MakeBranchInTreeS(treeS);
404     }
405     AliZDCMergedHit *MHit;
406     TClonesArray &sdigits = *fMergedHits;
407     Int_t imhit;
408     //Merged Hits loop
409     for(imhit=0; imhit<fNMergedhits; imhit++){
410        MHit = (AliZDCMergedHit*) fMergedHits->UncheckedAt(imhit);
411        new (sdigits[imhit]) AliZDCMergedHit(*MHit);
412     }
413     treeS->Fill();
414     treeS->AutoSave();
415   }
416   
417 }
418
419 //_____________________________________________________________________________
420 void AliZDC::SDigits2Digits()
421 {
422   if(!fMerger){ // Only digitization
423     printf("    ZDC digitization (no merging) \n");
424     fMerger = new AliZDCMerger();    
425     fMerger->Digitize(fNMergedhits, fMergedHits);
426
427     char hname[30];
428     sprintf(hname,"TreeD%d",gAlice->GetHeader()->GetEvent());
429     gAlice->TreeD()->Fill();
430     gAlice->TreeD()->AutoSave();
431     gAlice->TreeD()->Reset();  
432   }
433   else if(fMerger){     // Merging and digitization
434     printf("    ZDC merging and digitization\n");
435     fMerger->Digitize(fNMergedhits, fMergedHits);
436
437     // Digits tree
438
439     TTree *treeD = fLoader->TreeD();
440     if (treeD == 0x0)
441      {
442       Int_t retval = fLoader->LoadDigits();
443       if (retval)
444        {
445          Error("SDigits2Digits","Error while loading Digits");
446          return;
447        }
448       treeD = fLoader->TreeD();
449      }
450
451
452
453     if(!treeD){
454       printf("\n ERROR -> Can't find TreeD%d in background file\n",fMerger->EvNum());
455     }    
456     // Branch address
457     char branchDname[20];
458     sprintf(branchDname,"%s",GetName());
459     if(treeD && fDigits){
460       TBranch *branchD = treeD->GetBranch(branchDname);
461       if(branchD) branchD->SetAddress(&fDigits);
462       else if(!branchD) MakeBranchInTreeD(treeD);
463     }
464     treeD->Fill();
465     treeD->AutoSave();
466   }
467   
468   
469 }
470 //_____________________________________________________________________________
471 void AliZDC::Hits2Digits()
472 {
473     gAlice->Hits2SDigits();
474     gAlice->SDigits2Digits();
475 }
476
477 //_____________________________________________________________________________
478 void AliZDC::Digits2Reco()
479 {
480     printf("    Entering AliZDC::Digits2Reco\n");
481     AliDetector *ZDC  = gAlice->GetDetector("ZDC");
482     TClonesArray *ZDCdigits = ZDC->Digits();
483     
484     TTree *TD = fLoader->TreeD();
485     if (TD == 0x0)
486      {
487       Int_t retval = fLoader->LoadDigits();
488       if (retval)
489        {
490          Error("Digits2Reco","Error while loading Digits");
491          return;
492        }
493       TD = fLoader->TreeD();
494      }
495
496
497     if(TD){
498       char brname[20];
499       sprintf(brname,"%s",ZDC->GetName());
500       TBranch *br = TD->GetBranch(brname);
501       if(br) br->SetAddress(&ZDCdigits);
502     }
503     else if(!TD) printf("       ERROR -> TreeD NOT found in gAlice object\n");
504     
505     Int_t nt = (Int_t) (TD->GetEntries());
506     gAlice->ResetDigits();    
507     
508     AliZDCDigit *dig;
509     Int_t j, idig, ndigits, ZNraw=0, ZPraw=0, ZEMraw=0;
510     //  ---      Summing raw ADCs for each detector to obtain total light
511     for(j=0; j<nt; j++){
512       TD->GetEvent(j);
513       ndigits = ZDCdigits->GetEntries();
514       ZNraw=0;
515       ZPraw=0; 
516       ZEMraw=0;
517       //  ---  Loop over event digits
518       for(idig=0; idig<ndigits; idig++){
519          dig = (AliZDCDigit*) ZDCdigits->UncheckedAt(idig);
520          if(dig->GetSector(0) == 1)      ZNraw  += dig->GetADCValue();
521          else if(dig->GetSector(0) == 2) ZPraw  += dig->GetADCValue();
522          else if(dig->GetSector(0) == 3) ZEMraw += dig->GetADCValue();
523       } // Digits loop
524     } //  TreeD entries loop
525     printf("\n  ---     ZNraw = %d, ZPraw = %d, ZEMraw = %d\n",ZNraw, ZPraw, ZEMraw);
526     
527   //  ---      Pedestal subtraction
528   Int_t ZNcorr, ZPcorr, ZEMcorr, MeanPed=50;
529   ZNcorr  = ZNraw  - 5*MeanPed;
530   ZPcorr  = ZPraw  - 5*MeanPed;
531   ZEMcorr = ZEMraw - 2*MeanPed;
532   if(ZNcorr<0)  ZNcorr=0;
533   if(ZPcorr<0)  ZPcorr=0;
534   if(ZEMcorr<0) ZEMcorr=0;
535  printf("\n    ZNcorr = %d, ZPcorr = %d, ZEMcorr = %d\n",ZNcorr,ZPcorr,ZEMcorr);
536   
537   //  ---      ADCchannel -> photoelectrons
538   // NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7)
539   Float_t ZNphe, ZPphe, ZEMphe, ConvFactor = 0.064;
540   ZNphe  = ZNcorr/ConvFactor;
541   ZPphe  = ZPcorr/ConvFactor;
542   ZEMphe = ZEMcorr/ConvFactor;
543   printf("\n    ZNphe = %f, ZPphe = %f, ZEMphe = %f\n",ZNphe, ZPphe, ZEMphe);
544   
545   //  ---      Energy calibration
546   // Conversion factors for hadronic ZDCs goes from phe yield to TRUE incident 
547   //  energy (conversion from GeV to TeV is included); while for EM calos 
548   // conversion is from light yield to detected energy calculated by GEANT
549   // NB -> ZN and ZP conversion factors are constant since incident spectators
550   // have all the same energy, ZEM energy is obtained through a fit over the whole
551   // range of incident particle energies (obtained with full HIJING simulations) 
552   Float_t ZNenergy, ZPenergy, ZEMenergy, ZDCenergy;
553   Float_t ZNphexTeV=329., ZPphexTeV=369.;
554   ZNenergy  = ZNphe/ZNphexTeV;
555   ZPenergy  = ZPphe/ZPphexTeV;
556   ZDCenergy = ZNenergy+ZPenergy;
557   ZEMenergy = -4.81+0.3238*ZEMphe;
558   if(ZEMenergy<0) ZEMenergy=0;
559   printf("    ZNenergy = %f TeV, ZPenergy = %f TeV, ZDCenergy = %f GeV, "
560          "\n            ZEMenergy = %f TeV\n", ZNenergy, ZPenergy, 
561          ZDCenergy, ZEMenergy);
562   
563   if(ZDCenergy==0)
564     printf("\n\n        ###     ATTENZIONE!!! -> ev# %d: ZNenergy = %f TeV, ZPenergy = %f TeV, ZDCenergy = %f GeV, "
565          " ZEMenergy = %f TeV\n\n", fMerger->EvNum(), ZNenergy, ZPenergy, ZDCenergy, ZEMenergy); 
566   
567   //  ---      Number of incident spectator nucleons
568   Int_t NDetSpecN, NDetSpecP;
569   NDetSpecN = (Int_t) (ZNenergy/2.760);
570   NDetSpecP = (Int_t) (ZPenergy/2.760);
571   printf("\n    NDetSpecN = %d, NDetSpecP = %d\n",NDetSpecN, NDetSpecP);
572   
573   //  ---      Number of generated spectator nucleons and impact parameter
574   // --------------------------------------------------------------------------------------------------
575   // [1] ### Results in Chiara's PhD thesis -> 0<b<15 fm (Dec 2001)
576   /*// Fit results for neutrons (Nspectator n true vs. EZN)
577   TF1 *fZNCen = new TF1("fZNCen",
578       "(-2.116909+sqrt(2.116909*2.116909-4*(-0.00651)*(14.556798-x)))/(2*(-0.00651))",0.,158.5);
579   TF1 *fZNPer = new TF1("fZNPer",
580       "(-34.695134-sqrt(34.695134*34.695134-4*(-0.174780)*(-1562.283443-x)))/(2*(-0.174780))",0.,158.5);
581   // Fit results for protons (Nspectator p true vs. EZP)
582   TF1 *fZPCen = new TF1("fZPCen",
583       "(-1.3217+sqrt(1.3217*1.3217-4*(-0.007934)*(4.742873-x)))/(2*(-0.007934))",0.,58.91);
584   TF1 *fZPPer = new TF1("fZPPer",
585       "(-15.788267-sqrt(15.788267*15.788267-4*(-0.133359)*(-383.800673-x)))/(2*(-0.133359))",0.,58.91);
586   // Fit results for total number of spectators (Nspectators true vs. EZDC)
587   TF1 *fZDCCen = new TF1("fZDCCen",
588       "(-1.867335+sqrt(1.867335*1.867335-4*(-0.004119)*(19.100289-x)))/(2*(-0.004119))",0.,220.4);
589   TF1 *fZDCPer = new TF1("fZDCPer",
590       "(-22.429097-sqrt(22.429097*22.429097-4*(-0.072435)*(-1482.034526-x)))/(2*(-0.072435))",0.,220.4);*/
591   // --------------------------------------------------------------------------------------------------
592   // [1] ### Results from a new production  -> 0<b<18 fm (Apr 2002)
593   // Fit results for neutrons (Nspectator n true vs. EZN)
594   TF1 *fZNCen = new TF1("fZNCen",
595       "(-2.287920+sqrt(2.287920*2.287920-4*(-0.007629)*(11.921710-x)))/(2*(-0.007629))",0.,164.);
596   TF1 *fZNPer = new TF1("fZNPer",
597       "(-37.812280-sqrt(37.812280*37.812280-4*(-0.190932)*(-1709.249672-x)))/(2*(-0.190932))",0.,164.);
598   // Fit results for protons (Nspectator p true vs. EZP)
599   TF1 *fZPCen = new TF1("fZPCen",
600        "(-1.321353+sqrt(1.321353*1.321353-4*(-0.007283)*(3.550697-x)))/(2*(-0.007283))",0.,60.);
601   TF1 *fZPPer = new TF1("fZPPer",
602       "(-42.643308-sqrt(42.643308*42.643308-4*(-0.310786)*(-1402.945615-x)))/(2*(-0.310786))",0.,60.);
603   // Fit results for total number of spectators (Nspectators true vs. EZDC)
604   TF1 *fZDCCen = new TF1("fZDCCen",
605       "(-1.934991+sqrt(1.934991*1.934991-4*(-0.004080)*(15.111124-x)))/(2*(-0.004080))",0.,225.);
606   TF1 *fZDCPer = new TF1("fZDCPer",
607       "(-34.380639-sqrt(34.380639*34.380639-4*(-0.104251)*(-2612.189017-x)))/(2*(-0.104251))",0.,225.);
608   // --------------------------------------------------------------------------------------------------
609   // [1] ### Results in Chiara's PhD thesis -> 0<b<15 fm (Dec 2001)
610   /*// Fit results for b (b vs. EZDC)
611   //TF1 *fbCen = new TF1("fbCen","0.611543+0.052231*x-0.000112*x*x+0.000000374*x*x*x",0.,222.);
612   //TF1 *fbPer = new TF1("fbPer","16.552010-0.023866*x-0.00001*x*x",0.,222.);
613   TF1 *fbCen = new TF1("fbCen","0.612769+0.051929*x-0.0001074*x*x+0.0000003724*x*x*x",0.,225.);
614   TF1 *fbPer = new TF1("fbPer","16.6131016-0.026053*x+0.000006893*x*x",0.,225.);*/
615   // --------------------------------------------------------------------------------------------------
616   // [2] ### Results from a new production  -> 0<b<18 fm (Apr 2002)
617   TF1 *fbCen = new TF1("fbCen","-0.056923+0.079703*x-0.0004301*x*x+0.000001366*x*x*x",0.,220.);
618   TF1 *fbPer = new TF1("fbPer","17.943998-0.046846*x+0.000074*x*x",0.,220.);
619   // --------------------------------------------------------------------------------------------------
620   // Evaluating Nspectators and b from ZEM energy
621   // [1] ### Results in Chiara's PhD thesis -> 0<b<15 fm (Dec 2001)
622   /*TF1 *fZEMn  = new TF1("fZEMn","124.2-0.0566*x+0.000006014*x*x",0.,3500.);
623   TF1 *fZEMp  = new TF1("fZEMp","81.3-0.03834*x+0.000004359*x*x",0.,3500.);
624   TF1 *fZEMsp = new TF1("fZEMsp","205.6-0.09567*x+0.00001056*x*x",0.,3500.);
625   TF1 *fZEMb  = new TF1("fZEMb","15.8-0.02084*x+2.802e-5*x*x-2.007e-8*x*x*x+6.586e-12*x*x*x*x-8.042e-16*x*x*x*x*x",0.,3500.);*/
626   // --------------------------------------------------------------------------------------------------
627   // [2] ### Results from a new production  -> 0<b<18 fm (Apr 2002)
628   TF1 *fZEMn  = new TF1("fZEMn","126.2-0.05399*x+0.000005679*x*x",0.,4000.);
629   TF1 *fZEMp  = new TF1("fZEMp","82.49-0.03611*x+0.00000385*x*x",0.,4000.);
630   TF1 *fZEMsp = new TF1("fZEMsp","208.7-0.09006*x+0.000009526*x*x",0.,4000.);
631   TF1 *fZEMb  = new TF1("fZEMb","16.06-0.01633*x+1.44e-5*x*x-6.778e-9*x*x*x+1.438e-12*x*x*x*x-1.112e-16*x*x*x*x*x",0.,4000.);
632   
633   Int_t NGenSpecN=0, NGenSpecP=0, NGenSpec=0;
634   Double_t ImpPar=0;
635   // Cut value for Ezem (GeV)
636   // [1] ### Results in Chiara's PhD thesis -> 0<b<15 fm (Dec 2001)
637   //Float_t EZEMCut = 360.; 
638   // [2] ### Results from a new production  -> 0<b<18 fm (Apr 2002)
639   Float_t EZEMCut = 420.;
640   Float_t DeltaEZEMSup = 690.; 
641   Float_t DeltaEZEMInf = 270.; 
642   if(ZEMenergy > (EZEMCut+DeltaEZEMSup)){
643     NGenSpecN = (Int_t) (fZNCen->Eval(ZNenergy));
644     NGenSpecP = (Int_t) (fZPCen->Eval(ZPenergy));
645     NGenSpec  = (Int_t) (fZDCCen->Eval(ZDCenergy));
646     ImpPar    = fbCen->Eval(ZDCenergy);
647     //printf("    fZNCen = %f, fZPCen = %f, fZDCCen = %f\n",fZNCen->Eval(ZNenergy),
648     //            fZPCen->Eval(ZPenergy),fZDCCen->Eval(ZDCenergy));
649   }
650   else if(ZEMenergy < (EZEMCut-DeltaEZEMInf)){
651     NGenSpecN = (Int_t) (fZNPer->Eval(ZNenergy)); 
652     NGenSpecP = (Int_t) (fZPPer->Eval(ZPenergy));
653     NGenSpec  = (Int_t) (fZDCPer->Eval(ZDCenergy));
654     ImpPar    = fbPer->Eval(ZDCenergy);
655     //printf("    fZNPer = %f, fZPPer = %f, fZDCPer = %f\n",fZNPer->Eval(ZNenergy),
656     //            fZPPer->Eval(ZPenergy),fZDCPer->Eval(ZDCenergy));
657   }
658   else if(ZEMenergy >= (EZEMCut-DeltaEZEMInf) && ZEMenergy <= (EZEMCut+DeltaEZEMSup)){
659     NGenSpecN = (Int_t) (fZEMn->Eval(ZEMenergy));
660     NGenSpecP = (Int_t) (fZEMp->Eval(ZEMenergy));
661     NGenSpec  = (Int_t)(fZEMsp->Eval(ZEMenergy));
662     ImpPar    =  fZEMb->Eval(ZEMenergy);
663     //printf("    Nspec ZEM = %f, Nspec ZDC = %f\n",fZEMsp->Eval(ZNenergy),fZDCPer->Eval(ZDCenergy));
664   }
665   // [1] ### Results in Chiara's PhD thesis -> 0<b<15 fm (Dec 2001)
666   /*if(ZNenergy>158.5)  NGenSpecN = (Int_t) (fZEMn->Eval(ZEMenergy));
667   if(ZPenergy>58.91)  NGenSpecP = (Int_t) (fZEMp->Eval(ZEMenergy));
668   if(ZDCenergy>220.4) NGenSpec  = (Int_t)(fZEMsp->Eval(ZEMenergy));
669   if(ZDCenergy>225.)  ImpPar    =          fZEMb->Eval(ZEMenergy);*/
670   // [2] ### Results from a new production  -> 0<b<18 fm (Apr 2002)
671   if(ZNenergy>162.)  NGenSpecN = (Int_t) (fZEMn->Eval(ZEMenergy));
672   if(ZPenergy>59.75)  NGenSpecP = (Int_t) (fZEMp->Eval(ZEMenergy));
673   if(ZDCenergy>221.5) NGenSpec  = (Int_t)(fZEMsp->Eval(ZEMenergy));
674   if(ZDCenergy>220.)  ImpPar    =  fZEMb->Eval(ZEMenergy);
675   
676   if(NGenSpecN>125)    NGenSpecN=125;
677   else if(NGenSpecN<0) NGenSpecN=0;
678   if(NGenSpecP>82)     NGenSpecP=82;
679   else if(NGenSpecP<0) NGenSpecP=0;
680   if(NGenSpec>207)     NGenSpec=207;
681   else if(NGenSpec<0)  NGenSpec=0;
682   //printf("    NRecSpecN = %d, NRecSpecP = %d, NRecSpec = %d\n",NGenSpecN,NGenSpecP,NGenSpec);
683   
684   //  ---      Number of participants
685   Int_t NPart, NPartTot;
686   NPart = 207-NGenSpecN-NGenSpecP;
687   NPartTot = 207-NGenSpec;
688   //printf("    ###     NPart(ZP+ZN) = %d, NPart(ZDC) = %d, b = %f fm\n",NPart,NPartTot,ImpPar);
689   printf("      ###     NPart = %d, b = %f fm\n",NPartTot,ImpPar);
690   
691   //  ---     Writing RecPoints TCA
692   // Allocate the RecPoints TCA 
693   fRecPoints = new TClonesArray("AliZDCReco",1000);
694   AliZDCReco *reco = new AliZDCReco(ZNenergy,ZPenergy,ZDCenergy,ZEMenergy,
695               NDetSpecN,NDetSpecP,NGenSpecN,NGenSpecP,NGenSpec,NPartTot,ImpPar);
696   new((*fRecPoints)[fNRecPoints]) AliZDCReco(*reco);
697   //fNRecPoints++;
698   //fRecPoints->Dump();
699   delete reco;
700   
701   // TreeR
702   TTree *treeR = gAlice->TreeR();
703   if(!treeR) printf("\n ERROR -> Can't find TreeR%d in background file\n",fMerger->EvNum());
704   // Branch address
705   char branchRname[20];
706   sprintf(branchRname,"%s",GetName());
707   if(fRecPoints){
708     TBranch *branchR = treeR->GetBranch(branchRname);
709     if(branchR) branchR->SetAddress(&fRecPoints);
710     else if(!branchR) MakeBranchInTreeR(treeR);
711   }
712   treeR->Fill();
713   treeR->AutoSave();
714   treeR->Reset();
715 }
716
717