]>
Commit | Line | Data |
---|---|---|
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 | // Implementation version v1 of the EMCAL particle identifier | |
20 | ||
21 | //*-- Author: Yves Schutz (SUBATECH) | |
22 | // | |
23 | // --- ROOT system --- | |
24 | #include "TROOT.h" | |
25 | #include "TTree.h" | |
26 | #include "TFile.h" | |
27 | #include "TF2.h" | |
28 | #include "TCanvas.h" | |
29 | #include "TFolder.h" | |
30 | #include "TSystem.h" | |
31 | #include "TBenchmark.h" | |
32 | ||
33 | #include "TSystem.h" | |
34 | ||
35 | // --- Standard library --- | |
36 | ||
37 | #include <Riostream.h> | |
38 | ||
39 | // --- AliRoot header files --- | |
40 | ||
41 | #include "AliRun.h" | |
42 | #include "AliGenerator.h" | |
43 | #include "AliEMCAL.h" | |
44 | #include "AliEMCALPIDv1.h" | |
45 | #include "AliEMCALClusterizerv1.h" | |
46 | #include "AliEMCALTrackSegment.h" | |
47 | #include "AliEMCALTrackSegmentMakerv1.h" | |
48 | #include "AliEMCALRecParticle.h" | |
49 | #include "AliEMCALGeometry.h" | |
50 | #include "AliEMCALGetter.h" | |
51 | ||
52 | ClassImp( AliEMCALPIDv1) | |
53 | ||
54 | //____________________________________________________________________________ | |
55 | AliEMCALPIDv1::AliEMCALPIDv1():AliEMCALPID() | |
56 | { | |
57 | // default ctor | |
58 | ||
59 | InitParameters() ; | |
60 | fDefaultInit = kTRUE ; | |
61 | ||
62 | } | |
63 | ||
64 | //____________________________________________________________________________ | |
65 | AliEMCALPIDv1::AliEMCALPIDv1(const char * headerFile,const char * name, const Bool_t toSplit) | |
66 | :AliEMCALPID(headerFile, name,toSplit) | |
67 | { | |
68 | //ctor with the indication on where to look for the track segments | |
69 | ||
70 | InitParameters() ; | |
71 | ||
72 | Init() ; | |
73 | fDefaultInit = kFALSE ; | |
74 | ||
75 | } | |
76 | ||
77 | //____________________________________________________________________________ | |
78 | AliEMCALPIDv1::~AliEMCALPIDv1() | |
79 | { | |
80 | // dtor | |
81 | ||
82 | if (!fDefaultInit) { | |
83 | fSplitFile = 0 ; | |
84 | } | |
85 | } | |
86 | ||
87 | //____________________________________________________________________________ | |
88 | const TString AliEMCALPIDv1::BranchName() const | |
89 | { | |
90 | TString branchName(GetName() ) ; | |
91 | branchName.Remove(branchName.Index(Version())-1) ; | |
92 | return branchName ; | |
93 | } | |
94 | ||
95 | //____________________________________________________________________________ | |
96 | void AliEMCALPIDv1::Init() | |
97 | { | |
98 | // Make all memory allocations that are not possible in default constructor | |
99 | // Add the PID task to the list of EMCAL tasks | |
100 | ||
101 | if ( strcmp(GetTitle(), "") == 0 ) | |
102 | SetTitle("galice.root") ; | |
103 | ||
104 | TString branchname(GetName()) ; | |
105 | branchname.Remove(branchname.Index(Version())-1) ; | |
106 | AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(),branchname.Data(),fToSplit ) ; | |
107 | ||
108 | // gime->SetRecParticlesTitle(BranchName()) ; | |
109 | if ( gime == 0 ) { | |
110 | Error("Init", "Could not obtain the Getter object !" ) ; | |
111 | return ; | |
112 | } | |
113 | ||
114 | fSplitFile = 0 ; | |
115 | if(fToSplit){ | |
116 | //First - extract full path if necessary | |
117 | TString fileName(GetTitle()) ; | |
118 | Ssiz_t islash = fileName.Last('/') ; | |
119 | if(islash<fileName.Length()) | |
120 | fileName.Remove(islash+1,fileName.Length()) ; | |
121 | else | |
122 | fileName="" ; | |
123 | fileName+="EMCAL.RecData." ; | |
124 | if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){ | |
125 | fileName+=branchname.Data() ; | |
126 | fileName+="." ; | |
127 | } | |
128 | fileName+="root" ; | |
129 | fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data())); | |
130 | if(!fSplitFile) | |
131 | fSplitFile = TFile::Open(fileName.Data(),"update") ; | |
132 | } | |
133 | ||
134 | gime->PostPID(this) ; | |
135 | gime->PostRecParticles(branchname) ; | |
136 | ||
137 | } | |
138 | ||
139 | //____________________________________________________________________________ | |
140 | void AliEMCALPIDv1::InitParameters() | |
141 | { | |
142 | ||
143 | fRecParticlesInRun = 0 ; | |
144 | fNEvent = 0 ; | |
145 | fRecParticlesInRun = 0 ; | |
146 | TString pidName( GetName()) ; | |
147 | if (pidName.IsNull() ) | |
148 | pidName = "Default" ; | |
149 | pidName.Append(":") ; | |
150 | pidName.Append(Version()) ; | |
151 | SetName(pidName) ; | |
152 | fPi0Analysis = kFALSE ; | |
153 | } | |
154 | ||
155 | //____________________________________________________________________________ | |
156 | ||
157 | void AliEMCALPIDv1::Exec(Option_t * option) | |
158 | { | |
159 | //Steering method | |
160 | ||
161 | if( strcmp(GetName(), "")== 0 ) | |
162 | Init() ; | |
163 | ||
164 | if(strstr(option,"tim")) | |
165 | gBenchmark->Start("EMCALPID"); | |
166 | ||
167 | if(strstr(option,"print")) { | |
168 | Print("") ; | |
169 | return ; | |
170 | } | |
171 | AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; | |
172 | if(gime->BranchExists("RecParticles") ) | |
173 | return ; | |
174 | Int_t nevents = gime->MaxEvent() ; //(Int_t) gAlice->TreeE()->GetEntries() ; | |
175 | Int_t ievent ; | |
176 | ||
177 | ||
178 | for(ievent = 0; ievent < nevents; ievent++){ | |
179 | gime->Event(ievent,"R") ; | |
180 | ||
181 | MakeRecParticles() ; | |
182 | ||
183 | WriteRecParticles(ievent); | |
184 | ||
185 | if(strstr(option,"deb")) | |
186 | PrintRecParticles(option) ; | |
187 | ||
188 | //increment the total number of rec particles per run | |
189 | fRecParticlesInRun += gime->RecParticles(BranchName())->GetEntriesFast() ; | |
190 | ||
191 | } | |
192 | ||
193 | if(strstr(option,"tim")){ | |
194 | gBenchmark->Stop("EMCALPID"); | |
195 | Info("Exec", "took %f seconds for PID %f seconds per event", | |
196 | gBenchmark->GetCpuTime("EMCALPID"), | |
197 | gBenchmark->GetCpuTime("EMCALPID")/nevents) ; | |
198 | } | |
199 | } | |
200 | ||
201 | //____________________________________________________________________________ | |
202 | void AliEMCALPIDv1::MakeRecParticles(){ | |
203 | ||
204 | // Makes a RecParticle out of a TrackSegment | |
205 | ||
206 | AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; | |
207 | TObjArray * aECRecPoints = gime->ECALRecPoints() ; | |
208 | TObjArray * aPRRecPoints = gime->PRERecPoints() ; | |
209 | TObjArray * aHCRecPoints = gime->HCALRecPoints() ; | |
210 | TClonesArray * trackSegments = gime->TrackSegments() ; | |
211 | if ( !aECRecPoints || !aPRRecPoints || !aHCRecPoints || !trackSegments ) { | |
212 | Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ; | |
213 | } | |
214 | TClonesArray * recParticles = gime->RecParticles() ; | |
215 | recParticles->Clear(); | |
216 | ||
217 | TIter next(trackSegments) ; | |
218 | AliEMCALTrackSegment * ts ; | |
219 | Int_t index = 0 ; | |
220 | AliEMCALRecParticle * rp ; | |
221 | while ( (ts = (AliEMCALTrackSegment *)next()) ) { | |
222 | ||
223 | new( (*recParticles)[index] ) AliEMCALRecParticle() ; | |
224 | rp = (AliEMCALRecParticle *)recParticles->At(index) ; | |
225 | rp->SetTrackSegment(index) ; | |
226 | rp->SetIndexInList(index) ; | |
227 | ||
228 | AliEMCALTowerRecPoint * ecal = 0 ; | |
229 | if(ts->GetECIndex()>=0) | |
230 | ecal = dynamic_cast<AliEMCALTowerRecPoint *>(aECRecPoints->At(ts->GetECIndex())) ; | |
231 | ||
232 | AliEMCALTowerRecPoint * pre = 0 ; | |
233 | if(ts->GetPREIndex()>=0) | |
234 | pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRRecPoints->At(ts->GetPREIndex())) ; | |
235 | ||
236 | AliEMCALTowerRecPoint * hcal = 0 ; | |
237 | if(ts->GetHCIndex()>=0) | |
238 | hcal = dynamic_cast<AliEMCALTowerRecPoint *>(aHCRecPoints->At(ts->GetHCIndex())) ; | |
239 | ||
240 | // Now set type (reconstructed) of the particle | |
241 | ||
242 | // Choose the cluster energy range | |
243 | ||
244 | if (!ecal) { | |
245 | Fatal("MakeRecParticles", "-> emcal(%d) = %d", ts->GetECIndex(), ecal ) ; | |
246 | } | |
247 | ||
248 | Float_t e = ecal->GetEnergy() ; | |
249 | ||
250 | Float_t lambda[2] ; | |
251 | ecal->GetElipsAxis(lambda) ; | |
252 | ||
253 | if((lambda[0]>0.01) && (lambda[1]>0.01)){ | |
254 | // Looking PCA. Define and calculate the data (X), | |
255 | // introduce in the function X2P that gives the components (P). | |
256 | ||
257 | Float_t Spher = 0. ; | |
258 | Float_t Emaxdtotal = 0. ; | |
259 | ||
260 | if((lambda[0]+lambda[1])!=0) | |
261 | Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); | |
262 | ||
263 | Emaxdtotal=ecal->GetMaximalEnergy()/ecal->GetEnergy(); | |
264 | } | |
265 | ||
266 | // Float_t time = ecal->GetTime() ; | |
267 | ||
268 | //Set momentum, energy and other parameters | |
269 | Float_t encal = e; | |
270 | TVector3 dir(0., 0., 0.) ; | |
271 | dir.SetMag(encal) ; | |
272 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; | |
273 | rp->SetCalcMass(0); | |
274 | rp->Name(); //If photon sets the particle pdg name to gamma | |
275 | rp->SetProductionVertex(0,0,0,0); | |
276 | rp->SetFirstMother(-1); | |
277 | rp->SetLastMother(-1); | |
278 | rp->SetFirstDaughter(-1); | |
279 | rp->SetLastDaughter(-1); | |
280 | rp->SetPolarisation(0,0,0); | |
281 | index++ ; | |
282 | } | |
283 | ||
284 | } | |
285 | ||
286 | //____________________________________________________________________________ | |
287 | void AliEMCALPIDv1:: Print() | |
288 | { | |
289 | // Print the parameters used for the particle type identification | |
290 | ||
291 | TString message ; | |
292 | message = "\n=============== AliEMCALPID1 ================\n" ; | |
293 | message += "Making PID\n"; | |
294 | message += " Pricipal analysis file from 0.5 to 100 %s\n" ; | |
295 | message += " Name of parameters file %s\n" ; | |
296 | message += " Matrix of Parameters: 9x4\n" ; | |
297 | message += " RCPV 2x3 rows x and z, columns function cut parameters\n" ; | |
298 | message += " TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n" ; | |
299 | message += " PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n" ; | |
300 | message += " Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n" ; | |
301 | Info("Print", message.Data(), fFileName.Data(), fFileNamePar.Data() ) ; | |
302 | } | |
303 | ||
304 | //____________________________________________________________________________ | |
305 | void AliEMCALPIDv1::WriteRecParticles(Int_t event) | |
306 | { | |
307 | ||
308 | AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; | |
309 | ||
310 | TClonesArray * recParticles = gime->RecParticles() ; | |
311 | recParticles->Expand(recParticles->GetEntriesFast() ) ; | |
312 | TTree * treeR ; | |
313 | ||
314 | if(fToSplit){ | |
315 | if(!fSplitFile) | |
316 | return ; | |
317 | fSplitFile->cd() ; | |
318 | char name[10] ; | |
319 | sprintf(name,"%s%d", "TreeR",event) ; | |
320 | treeR = dynamic_cast<TTree*>(fSplitFile->Get(name)); | |
321 | } | |
322 | else{ | |
323 | treeR = gAlice->TreeR(); | |
324 | } | |
325 | ||
326 | if(!treeR){ | |
327 | gAlice->MakeTree("R", fSplitFile); | |
328 | treeR = gAlice->TreeR() ; | |
329 | } | |
330 | ||
331 | //First rp | |
332 | Int_t bufferSize = 32000 ; | |
333 | TBranch * rpBranch = treeR->Branch("EMCALRP",&recParticles,bufferSize); | |
334 | rpBranch->SetTitle(BranchName()); | |
335 | ||
336 | ||
337 | //second, pid | |
338 | Int_t splitlevel = 0 ; | |
339 | AliEMCALPIDv1 * pid = this ; | |
340 | TBranch * pidBranch = treeR->Branch("AliEMCALPID","AliEMCALPIDv1",&pid,bufferSize,splitlevel); | |
341 | pidBranch->SetTitle(BranchName()); | |
342 | ||
343 | rpBranch->Fill() ; | |
344 | pidBranch->Fill() ; | |
345 | ||
346 | treeR->AutoSave() ; //Write(0,kOverwrite) ; | |
347 | if(gAlice->TreeR()!=treeR){ | |
348 | treeR->Delete(); | |
349 | } | |
350 | } | |
351 | ||
352 | //____________________________________________________________________________ | |
353 | void AliEMCALPIDv1::PrintRecParticles(Option_t * option) | |
354 | { | |
355 | // Print table of reconstructed particles | |
356 | ||
357 | AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; | |
358 | ||
359 | TClonesArray * recParticles = gime->RecParticles(BranchName()) ; | |
360 | ||
361 | TString message ; | |
362 | message = "\nevent " ; | |
363 | message += gAlice->GetEvNumber() ; | |
364 | message += " found " ; | |
365 | message += recParticles->GetEntriesFast(); | |
366 | message += " RecParticles\n" ; | |
367 | ||
368 | if(strstr(option,"all")) { // printing found TS | |
369 | message += "\n PARTICLE Index \n" ; | |
370 | ||
371 | Int_t index ; | |
372 | for (index = 0 ; index < recParticles->GetEntries() ; index++) { | |
373 | AliEMCALRecParticle * rp = (AliEMCALRecParticle * ) recParticles->At(index) ; | |
374 | message += "\n" ; | |
375 | message += rp->Name().Data() ; | |
376 | message += " " ; | |
377 | message += rp->GetIndexInList() ; | |
378 | message += " " ; | |
379 | message += rp->GetType() ; | |
380 | } | |
381 | } | |
382 | Info("Print", message.Data() ) ; | |
383 | } | |
384 | ||
385 | ||
386 |