Additional protection
[u/mrichter/AliRoot.git] / FMD / AliFMDReconstructor.cxx
CommitLineData
37c55dc0 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//_________________________________________________________________________
121a60bd 17// This is a class that constructs ReconstParticles (reconstructed particles)
37c55dc0 18// out of Digits
19//
20//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
21//////////////////////////////////////////////////////////////////////////////
22
dc8af42e 23// --- ROOT system ---
24#include "TTask.h"
25#include "TTree.h"
26#include "TSystem.h"
27#include "TFile.h"
37c55dc0 28#include "TROOT.h"
29#include "TFolder.h"
46501dfb 30#include "TH2F.h"
37c55dc0 31
dc8af42e 32// --- Standard library ---
37c55dc0 33#include <stdlib.h>
88cb7938 34#include <Riostream.h>
dc8af42e 35
36// --- AliRoot header files ---
37
88cb7938 38#include "AliRunLoader.h"
39#include "AliLoader.h"
40
dc8af42e 41#include "AliFMDdigit.h"
1a42d4f2 42#include "AliFMDhit.h"
dc8af42e 43#include "AliFMDReconstParticles.h"
44#include "AliFMD.h"
45#include "AliFMDv1.h"
121a60bd 46#include "AliFMDReconstructor.h"
dc8af42e 47#include "AliRun.h"
88cb7938 48#include "AliConfig.h"
1a42d4f2 49#include "AliHeader.h"
50#include "AliGenEventHeader.h"
88cb7938 51
121a60bd 52ClassImp(AliFMDReconstructor)
dc8af42e 53
54
88cb7938 55//____________________________________________________________________________
dc8af42e 56
121a60bd 57void AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
dc8af42e 58{
88cb7938 59 //Collects all digits in the same active volume into number of particles
a39f3926 60 /*
88cb7938 61 Reconstruct number of particles
a39f3926 62 in given group of pads for given FMDvolume
63 determine by numberOfVolume ,
64 numberOfMinSector,numberOfMaxSector,
65 numberOfMinRing, numberOgMaxRing
121a60bd 66 Reconstructor method choose dependence on number of empty pads
a39f3926 67 */
68
69
3d44ce66 70#ifdef DEBUG
71 Info("Exec ","Start");
72#endif
73
a39f3926 74 Int_t const knumVolumes=5;
46501dfb 75 Int_t padADC;
1a42d4f2 76 Float_t eta, etain,etaout,rad,theta;
cb1df35e 77 Int_t ivol, iSector, iRing;
1a42d4f2 78 Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
79 Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
2ad713c4 80 Float_t z[5]={-62.8, -75.2, 83.4, 75.2, 340.};
46501dfb 81 Int_t numberOfRings[5]={512,256,512,256,512};
a39f3926 82 Int_t numberOfSectors[5]= {20,40,20,40,20};
83 Int_t numberOfEtaIntervals[5];
cb1df35e 84 // number of ring for boundary 0.1 eta
cb1df35e 85
88cb7938 86
121a60bd 87 if (runLoader == 0x0)
88cb7938 88 {
89 Error("Exec","Run Loader loader is NULL - Session not opened");
90 return;
91 }
3d44ce66 92
121a60bd 93 AliLoader* plFMD = runLoader->GetLoader("FMDLoader");
3d44ce66 94 if (plFMD == 0x0)
88cb7938 95 {
121a60bd 96 Fatal("AliFMDReconstructor","Can not find FMD (loader) in specified event");
88cb7938 97 return;//never reached
98 }
99
121a60bd 100 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
101 if (!runLoader->TreeE()) runLoader->LoadHeader();
3d44ce66 102
46501dfb 103 TDirectory* cwd = gDirectory;
104 gDirectory = 0x0;
3d44ce66 105 Text_t buf1[20];
46501dfb 106 TH2F* hTotal[10];
107 for (Int_t j=1; j<=5; j++){
108 sprintf(buf1,"hTotal%d",j);
109
110 hTotal[j] = new TH2F(buf1," Number of primary particles ",
111 numberOfSectors[j-1],1,numberOfSectors[j-1],
112 numberOfRings[j-1],1,numberOfRings[j-1]);
113 }
114 gDirectory = cwd;
115
121a60bd 116 plFMD->LoadRecPoints("RECREATE");
3d44ce66 117 Int_t retval=0;
121a60bd 118 Int_t nevents=Int_t (runLoader->TreeE()->GetEntries());
3d44ce66 119#ifdef DEBUG
121a60bd 120 cout<<" nevents "<<nevents<<endl;
3d44ce66 121#endif
121a60bd 122 for(Int_t ievent=0;ievent<nevents;ievent++)
dc8af42e 123 {
3d44ce66 124#ifdef DEBUG
125 cout<<" *** event "<<ievent<<endl;
126#endif
121a60bd 127 runLoader->GetEvent(ievent) ;
3d44ce66 128 //event z-vertex for correction eta-rad dependence
121a60bd 129 AliHeader *header = runLoader->GetHeader();
3d44ce66 130 AliGenEventHeader* genHeader = header->GenEventHeader();
131 TArrayF *o = new TArrayF(3);
27301b0a 132 if (genHeader) genHeader->PrimaryVertex(*o);
3d44ce66 133 Float_t zVertex=o->At(2);
88cb7938 134
3d44ce66 135 for (Int_t i=0; i<5; i++)
46501dfb 136 hTotal[i+1]->Reset();
137
3d44ce66 138 retval = plFMD->LoadDigits("READ");
46501dfb 139 if (retval)
140 {
141 Error("Exec","Error occured while loading digits. Exiting.");
142 return;
143 }
144
3d44ce66 145 AliFMD * fFMD = (AliFMD *)gAlice->GetDetector("FMD");
46501dfb 146 TClonesArray *fReconParticles=fFMD->ReconParticles();
147 TClonesArray *fDigits=fFMD->Digits();
1a42d4f2 148
3d44ce66 149 TTree* treeD = plFMD->TreeD();
46501dfb 150 if (treeD == 0x0)
151 {
152 Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
153 return;
154 }
155
156 TBranch *brDigits=0;
157
158 brDigits=treeD->GetBranch("FMD");
3d44ce66 159
46501dfb 160 if (brDigits) {
161 brDigits->SetAddress(&fDigits);
46501dfb 162 }else{
163 cerr<<"EXEC Branch FMD digits not found"<<endl;
164 return;
165 }
166
3d44ce66 167 if(plFMD->TreeR()==0) plFMD->MakeTree("R");
88cb7938 168
dc8af42e 169 //Make branches
a39f3926 170 fFMD->MakeBranch("R");
1a42d4f2 171
37c55dc0 172
46501dfb 173 Int_t zeroADC=6;
3d44ce66 174 // read Digits
1a42d4f2 175 AliFMDdigit *fmdDigit;
a39f3926 176 if (fFMD)
dc8af42e 177 {
3d44ce66 178 plFMD->TreeD()->GetEvent(0);
46501dfb 179
180 Int_t nDigits=fDigits->GetEntries();
46501dfb 181 Int_t recParticles[6];
182 Int_t nRecPart=0 ;
3d44ce66 183 Int_t zeroPads=0;
184 Int_t numberOfPads=0;
185 Int_t pedestal;
186 Float_t channelWidth=(22400*50)/1024;
187 for (Int_t digit=0;digit<nDigits;digit++)
188 {
189 fmdDigit=(AliFMDdigit*)fDigits->UncheckedAt(digit);
190 ivol=fmdDigit->Volume();
191 iSector=fmdDigit->NumberOfSector();
192 iRing=fmdDigit->NumberOfRing();
193 pedestal=Int_t(gRandom->Gaus(500,250));
194 padADC= fmdDigit->ADCsignal()
195 -Int_t(Float_t(pedestal)/channelWidth);
196 if (padADC<0) padADC=0;
46501dfb 197 hTotal[ivol]->Fill(iSector,iRing,padADC);
3d44ce66 198 } //digit loop
199
200 //reconstruct multiplicity in 0.1 eta according Poisson distribution
46501dfb 201
3d44ce66 202 Int_t rmin=0; Int_t rmax=0;
203 Int_t smin=0; Int_t smax=0;
204 for (ivol=0; ivol<knumVolumes; ivol++)
205 {
206 Float_t ring2number=Float_t (numberOfRings[ivol])/
207 (rout[ivol]-rin[ivol]);
208 Float_t realZ=zVertex+z[ivol];
209 theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
210 etain = - TMath::Log( TMath::Tan(theta/2.));
211 theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
212 etaout=- TMath::Log( TMath::Tan(theta/2.));
213 numberOfEtaIntervals[ivol]=Int_t((etaout-etain)*10)-1;
214 eta=etain;
215 for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++)
216 {
217 theta = 2.*TMath::ATan (TMath::Exp (-eta));
218 Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
219 rmax= Int_t ( (radmin-rin[ivol])*ring2number+0.5);
220 eta=eta+0.1;
221 theta = 2.*TMath::ATan (TMath::Exp (-eta));
222 rad = TMath::Abs(realZ) * (TMath::Tan (theta));
223 rmin=Int_t( (rad-rin[ivol])*ring2number+0.5);
224
225 zeroPads=0;
226 smin=0;
227 smax=numberOfSectors[ivol];
228 numberOfPads=(rmax-rmin)*(smax-smin);
229 for (Int_t iring=rmin; iring<rmax; iring++)
230 {
231 for
232 (Int_t isector=0;
233 isector<numberOfSectors[ivol];
234 isector++)
235 {
236 if(hTotal[ivol+1]->GetBinContent(isector+1,iring+1)
237 <zeroADC) zeroPads++;}
238
239 } //ring
240 Float_t numberOfPads=Float_t(smax-smin)*Float_t(rmax-rmin);
241
242 Double_t lambda=-TMath::Log(Double_t(zeroPads)/numberOfPads);
243 Int_t fRecon=Int_t (lambda*numberOfPads+0.5);
244 recParticles[0]=ivol+1;
245 recParticles[1]=smin;
246 recParticles[2]=smax;
247 recParticles[3]=rmin;
248 recParticles[4]=rmax;
249 recParticles[5]= fRecon;
250 new((*fReconParticles)[nRecPart++]) AliFMDReconstParticles(recParticles);
251
46501dfb 252
1a42d4f2 253 } // eta
254 } // volume
cb1df35e 255
cb1df35e 256 }//if FMD
3d44ce66 257 plFMD->TreeR()->Reset();
258 plFMD->TreeR()->Fill();
259 plFMD->WriteRecPoints("OVERWRITE");
260 plFMD->UnloadDigits();
dc8af42e 261 } //event loop
121a60bd 262 plFMD->UnloadRecPoints();
3d44ce66 263#ifdef DEBUG
264 Info(" Exec"," finished");
265#endif
266 // delete hTotal[10];
267
88cb7938 268}
dc8af42e 269
121a60bd 270
271//_____________________________________________________________________________
272void AliFMDReconstructor::FillESD(AliRunLoader* /*runLoader*/,
273 AliESD* /*esd*/) const
274{
275// nothing to be done
276
277}
278