1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //_________________________________________________________________________
17 // This is a TTask that constructs ReconstParticles (reconstructed particles)
20 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
21 //////////////////////////////////////////////////////////////////////////////
23 // --- ROOT system ---
32 // --- Standard library ---
34 #include <Riostream.h>
36 // --- AliRoot header files ---
38 #include "AliRunLoader.h"
39 #include "AliLoader.h"
41 #include "AliFMDdigit.h"
42 #include "AliFMDhit.h"
43 #include "AliFMDReconstParticles.h"
46 #include "AliFMDReconstruction.h"
48 #include "AliConfig.h"
49 #include "AliHeader.h"
50 #include "AliGenEventHeader.h"
52 ClassImp(AliFMDReconstruction)
55 //____________________________________________________________________________
57 AliFMDReconstruction::AliFMDReconstruction():TTask("AliFMDReconstruction","")
59 fNevents = 0 ; // Number of events to rreconnstraction, 0 means all events in current file
63 //____________________________________________________________________________
65 AliFMDReconstruction::AliFMDReconstruction(AliRunLoader* rl):TTask("AliFMDReconstruction","")
70 Fatal("AliFMDReconstruction","Argument AliRunLoader* is null!");
74 fNevents = 0 ; // Number of events to rreconnstraction, 0 means all events in current file
77 AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
80 Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
81 return;//never reached
83 //add Task to //root/Tasks folder
84 gime->PostReconstructioner(this);
87 //____________________________________________________________________________
89 AliFMDReconstruction::~AliFMDReconstruction()
93 //____________________________________________________________________________
95 void AliFMDReconstruction::Exec()
97 //Collects all digits in the same active volume into number of particles
99 Reconstruct number of particles
100 in given group of pads for given FMDvolume
101 determine by numberOfVolume ,
102 numberOfMinSector,numberOfMaxSector,
103 numberOfMinRing, numberOgMaxRing
104 Reconstruction method choose dependence on number of empty pads
108 cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
109 Int_t const knumVolumes=5;
111 Float_t eta, etain,etaout,rad,theta;
112 Int_t ivol, iSector, iRing;
113 Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
114 Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
115 Float_t z[5]={-62.8, -75.2, 83.4, 75.2, 340.};
116 Int_t numberOfRings[5]={512,256,512,256,512};
117 Int_t numberOfSectors[5]= {20,40,20,40,20};
118 Int_t numberOfEtaIntervals[5];
119 // number of ring for boundary 0.1 eta
122 if (fRunLoader == 0x0)
124 Error("Exec","Run Loader loader is NULL - Session not opened");
129 AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
132 Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
133 return;//never reached
136 fRunLoader->LoadgAlice();
137 fRunLoader->LoadHeader();
139 TDirectory* cwd = gDirectory;
143 for (Int_t j=1; j<=5; j++){
144 sprintf(buf1,"hTotal%d",j);
146 hTotal[j] = new TH2F(buf1," Number of primary particles ",
147 numberOfSectors[j-1],1,numberOfSectors[j-1],
148 numberOfRings[j-1],1,numberOfRings[j-1]);
153 if(fNevents == 0) fNevents=fRunLoader->TreeE()->GetEntries();
154 cout<<" fNevents "<<fNevents<<endl;
155 for(Int_t ievent=0;ievent<fNevents;ievent++)
157 fRunLoader->GetEvent(ievent) ;
159 cout<<" ievent "<<ievent<<endl;
160 for (Int_t i=0; i<5; i++)
161 hTotal[i+1]->Reset();
163 retval = gime->LoadDigits("READ");
166 Error("Exec","Error occured while loading digits. Exiting.");
170 AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
171 TClonesArray *fReconParticles=fFMD->ReconParticles();
172 TClonesArray *fDigits=fFMD->Digits();
174 TTree* treeD = gime->TreeD();
177 Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
183 brDigits=treeD->GetBranch("FMD");
184 cout<<" brDigits "<<brDigits<<endl;
186 brDigits->SetAddress(&fDigits);
187 // fFMD->SetHitsAddressBranch(brHits);
189 cerr<<"EXEC Branch FMD digits not found"<<endl;
193 if(gime->TreeR()==0) gime->MakeTree("R");
196 fFMD->MakeBranch("R");
201 AliFMDdigit *fmdDigit;
204 gime->TreeD()->GetEvent(0);
206 Int_t nDigits=fDigits->GetEntries();
207 cout<<" nDigits "<<nDigits<<endl;
208 Int_t recParticles[6];
211 Int_t numberOfPads=0; //To avoid warning
213 Float_t channelWidth=(22400*50)/1024;
214 for (Int_t digit=0;digit<nDigits;digit++)
216 fmdDigit=(AliFMDdigit*)fDigits->UncheckedAt(digit);
217 ivol=fmdDigit->Volume();
218 iSector=fmdDigit->NumberOfSector();
219 iRing=fmdDigit->NumberOfRing();
220 pedestal=Int_t(gRandom->Gaus(500,250));
221 padADC= fmdDigit->ADCsignal()
222 -Int_t(Float_t(pedestal)/channelWidth);
223 if (padADC<0) padADC=0;
224 hTotal[ivol]->Fill(iSector,iRing,padADC);
226 Int_t rmin=0; Int_t rmax=0; //To avoid warning
227 Int_t smin=0; Int_t smax=0; //To avoid warning
228 AliHeader *header = fRunLoader->GetHeader();
229 AliGenEventHeader* genHeader = header->GenEventHeader();
230 TArrayF *o = new TArrayF(3);
231 genHeader->PrimaryVertex(*o);
232 Float_t zVertex=o->At(2);
233 for (ivol=0; ivol<knumVolumes; ivol++)
235 Float_t ring2number=Float_t (numberOfRings[ivol])/
236 (rout[ivol]-rin[ivol]);
237 Float_t realZ=zVertex+z[ivol];
238 theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
239 etain = - TMath::Log( TMath::Tan(theta/2.));
240 theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
241 etaout=- TMath::Log( TMath::Tan(theta/2.));
242 numberOfEtaIntervals[ivol]=(etaout-etain)*10-1;
244 for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++)
246 theta = 2.*TMath::ATan (TMath::Exp (-eta));
247 Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
248 rmax= Int_t ( (radmin-rin[ivol])*ring2number+0.5);
249 cout<<ivol<<" "<<" eta "<<eta<<" radmin "<<radmin<<
251 " Rmax "<<rmax<<" "<<endl;;
253 theta = 2.*TMath::ATan (TMath::Exp (-eta));
254 rad = TMath::Abs(realZ) * (TMath::Tan (theta));
255 rmin=Int_t( (rad-rin[ivol])*ring2number+0.5);
256 cout<<"eta "<<eta<<" rad "<<rad<<" Rmin "<<rmin<<endl;
257 // if(ivol==2&&e1==13) rmin=0;
260 smax=numberOfSectors[ivol];
261 numberOfPads=(rmax-rmin)*(smax-smin);
262 for (Int_t iring=rmin; iring<rmax; iring++)
266 isector<numberOfSectors[ivol];
269 if(hTotal[ivol+1]->GetBinContent(isector+1,iring+1)
270 <zeroADC) zeroPads++;}
273 Float_t numberOfPads=Float_t(smax-smin)*Float_t(rmax-rmin);
275 cout<<" zero "<<zeroPads++<<" pads "<<numberOfPads;
276 Double_t lambda=-TMath::Log(Double_t(zeroPads)/numberOfPads);
277 Int_t fRecon=(lambda*numberOfPads+0.5);
280 (Float_t)zeroPads/(Float_t)numberOfPads;
281 cout<<" zerosRatio "<<zerosRatio<<" recon "<<fRecon<<endl;
282 recParticles[0]=ivol+1;
283 recParticles[1]=smin;
284 recParticles[2]=smax;
285 recParticles[3]=rmin;
286 recParticles[4]=rmax;
287 recParticles[5]= fRecon;
288 new((*fReconParticles)[nRecPart++])
289 AliFMDReconstParticles(recParticles);
296 gime->TreeR()->Reset();
297 gime->TreeR()->Fill();
298 gime->WriteRecPoints("OVERWRITE");
300 cout<<"\nAliFMDReconstruction::Exec finished"<<endl;