6ad0bfa0 |
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 | |
b2a60966 |
16 | /* $Id$ */ |
17 | |
6ad0bfa0 |
18 | //_________________________________________________________________________ |
b2a60966 |
19 | // Implementation version v1 of the PHOS particle identifier |
7acf6008 |
20 | // Particle identification based on the |
21 | // - CPV information, |
a4e98857 |
22 | // - Preshower information (in MIXT or GPS2 geometries) |
7acf6008 |
23 | // - shower width. |
a4e98857 |
24 | // |
25 | // CPV or Preshower clusters should be closer in PHOS plane than fCpvEmcDistance (in cm). |
26 | // This parameter can be set by method SetCpvtoEmcDistanceCut(Float_t cut) |
7acf6008 |
27 | // |
28 | // One can set desirable ID method by the function SetIdentificationMethod(option). |
a4e98857 |
29 | // Presently the following options can be used together or separately : |
7acf6008 |
30 | // - "disp": use dispersion cut on shower width |
31 | // (width can be set by method SetDispersionCut(Float_t cut) |
a4e98857 |
32 | // - "ell" : use cut on the axis of the ellipse, drawn around shower |
7acf6008 |
33 | // (this cut can be changed by SetShowerProfileCut(char* formula), |
34 | // where formula - any function of two variables f(lambda[0],lambda[1]). |
35 | // Shower is considered as EM if f() > 0 ) |
36 | // One can visualize current cuts calling method PlotDispersionCuts(). |
37 | // |
a4e98857 |
38 | // use case: |
39 | // root [0] AliPHOSPIDv1 * p1 = new AliPHOSPIDv1("galice.root") |
40 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
41 | // root [1] p1->SetIdentificationMethod("disp ellipse") |
42 | // root [2] p1->ExecuteTask() |
43 | // root [3] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","ts1") |
44 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
7acf6008 |
45 | // // reading headers from file galice1.root and TrackSegments |
46 | // // with title "ts1" |
a4e98857 |
47 | // root [4] p2->SetRecParticlesBranch("rp1") |
7acf6008 |
48 | // // set file name for the branch RecParticles |
a4e98857 |
49 | // root [5] p2->ExecuteTask("deb all time") |
7acf6008 |
50 | // // available options |
51 | // // "deb" - prints # of reconstructed particles |
52 | // // "deb all" - prints # and list of RecParticles |
53 | // // "time" - prints benchmarking results |
54 | // |
55 | //*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) & |
56 | // Dmitri Peressounko (SUBATECH & Kurchatov Institute) |
57 | // Complitely redesined by Dmitri Peressounko, March 2001 |
6ad0bfa0 |
58 | |
59 | // --- ROOT system --- |
7acf6008 |
60 | #include "TROOT.h" |
61 | #include "TTree.h" |
62 | #include "TFile.h" |
63 | #include "TF2.h" |
64 | #include "TFormula.h" |
65 | #include "TCanvas.h" |
66 | #include "TFolder.h" |
67 | #include "TSystem.h" |
68 | #include "TBenchmark.h" |
6ad0bfa0 |
69 | // --- Standard library --- |
70 | |
de9ec31b |
71 | #include <iostream.h> |
7acf6008 |
72 | #include <iomanip.h> |
6ad0bfa0 |
73 | |
74 | // --- AliRoot header files --- |
75 | |
7acf6008 |
76 | #include "AliRun.h" |
77 | #include "AliGenerator.h" |
78 | #include "AliPHOS.h" |
26d4b141 |
79 | #include "AliPHOSPIDv1.h" |
7acf6008 |
80 | #include "AliPHOSClusterizerv1.h" |
6ad0bfa0 |
81 | #include "AliPHOSTrackSegment.h" |
7acf6008 |
82 | #include "AliPHOSTrackSegmentMakerv1.h" |
6ad0bfa0 |
83 | #include "AliPHOSRecParticle.h" |
84 | |
26d4b141 |
85 | ClassImp( AliPHOSPIDv1) |
6ad0bfa0 |
86 | |
1cb7c1ee |
87 | //____________________________________________________________________________ |
88 | AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID() |
89 | { |
a4e98857 |
90 | // default ctor |
7acf6008 |
91 | fIsInitialized = kFALSE ; |
92 | } |
93 | |
94 | //____________________________________________________________________________ |
95 | AliPHOSPIDv1::AliPHOSPIDv1(const char * headeFile,const char * tsBranchTitle):AliPHOSPID() |
96 | { |
a4e98857 |
97 | //ctor with the indication on where to look for the track segments |
dd5c4038 |
98 | |
7acf6008 |
99 | fHeaderFileName = headeFile ; |
100 | |
101 | fTSTitle = tsBranchTitle ; |
102 | |
103 | SetName("AliPHOSPID") ; |
104 | SetTitle("version1") ; |
2bd5457f |
105 | fIsInitialized = kFALSE ; |
7acf6008 |
106 | |
2bd5457f |
107 | Init() ; |
7acf6008 |
108 | |
109 | } |
110 | //____________________________________________________________________________ |
111 | AliPHOSPIDv1::~AliPHOSPIDv1() |
112 | { |
a4e98857 |
113 | //dtor |
7acf6008 |
114 | } |
115 | //____________________________________________________________________________ |
116 | void AliPHOSPIDv1::Init() |
117 | { |
a4e98857 |
118 | // Make all memory allocations that are not possible in default constructor |
2bd5457f |
119 | // Add the PID task to the list of PHOS tasks |
120 | |
7acf6008 |
121 | if(!fIsInitialized){ |
122 | if(fHeaderFileName.IsNull()) |
123 | fHeaderFileName = "galice.root" ; |
124 | |
125 | TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ; |
126 | |
127 | if(file == 0){ |
2bd5457f |
128 | if(fHeaderFileName.Contains("rfio")) // if we read file using HPSS |
129 | file = TFile::Open(fHeaderFileName.Data(),"update") ; |
130 | else |
131 | file = new TFile(fHeaderFileName.Data(),"update") ; |
7acf6008 |
132 | gAlice = (AliRun *) file->Get("gAlice") ; |
133 | } |
134 | |
135 | AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ; |
136 | fGeom = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() ); |
137 | |
138 | fTrackSegments = new TClonesArray("AliPHOSTrackSegment",1) ; |
139 | fTSMaker = new AliPHOSTrackSegmentMakerv1() ; |
140 | fEmcRecPoints = new TObjArray(1) ; |
141 | fCpvRecPoints = new TObjArray(1) ; |
142 | fClusterizer = new AliPHOSClusterizerv1() ; |
143 | fRecParticles = new TClonesArray("AliPHOSRecParticle",100) ; |
144 | |
145 | fFormula = new TFormula("LambdaCuts","(x>1)*(x<2.5)*(y>0)*(y<x)") ; |
146 | |
2bd5457f |
147 | fDispersion = 2.0 ; |
7acf6008 |
148 | fCpvEmcDistance = 3.0 ; |
149 | |
2bd5457f |
150 | //add Task to //YSAlice/tasks/Reconstructioner/PHOS |
151 | TFolder * alice = (TFolder*)gROOT->GetListOfBrowsables()->FindObject("YSAlice") ; |
152 | TTask * aliceRe = (TTask*)alice->FindObject("tasks/Reconstructioner") ; |
153 | TTask * phosRe = (TTask*)aliceRe->GetListOfTasks()->FindObject("PHOS") ; |
154 | phosRe->Add(this) ; |
155 | |
7acf6008 |
156 | fIsInitialized = kTRUE ; |
157 | } |
1cb7c1ee |
158 | } |
7acf6008 |
159 | //____________________________________________________________________________ |
160 | Bool_t AliPHOSPIDv1::ReadTrackSegments() |
161 | { |
f035f6ce |
162 | // Reads TrackSegments an extracts the title of the RecPoints |
a4e98857 |
163 | // branch from which TS were made of. |
f035f6ce |
164 | // Then reads both TrackSegments and RecPoints. |
7acf6008 |
165 | |
f035f6ce |
166 | //Fist read Track Segment Branch and extract RecPointsBranch from fTSMaker |
7acf6008 |
167 | fTrackSegments->Clear() ; |
168 | fEmcRecPoints->Clear() ; |
169 | fCpvRecPoints->Clear() ; |
170 | fRecParticles->Clear() ; |
171 | |
172 | gAlice->GetEvent(fNEvent) ; |
173 | |
174 | TTree * treeR = gAlice->TreeR() ; |
175 | |
176 | if(treeR==0){ |
177 | char treeName[20]; |
178 | sprintf(treeName,"TreeR%d",fNEvent); |
179 | cout << "Error in AliPHOSClusterizerv1 : no "<<treeName << endl ; |
180 | cout << " Do nothing " << endl ; |
181 | return kFALSE ; |
182 | } |
183 | |
184 | //first read TSMaker branch and extract information about RecPoints Branches |
185 | TBranch * tsMakerBranch = 0; |
186 | TBranch * tsBranch = 0; |
187 | |
188 | TObjArray * branches = treeR->GetListOfBranches() ; |
189 | Int_t ibranch; |
190 | Bool_t tsMakerNotFound = kTRUE ; |
191 | Bool_t tsNotFound = kTRUE ; |
192 | |
193 | for(ibranch = 0;(ibranch <branches->GetEntries())&&(tsMakerNotFound||tsNotFound);ibranch++){ |
194 | if(tsMakerNotFound){ |
195 | tsMakerBranch=(TBranch *) branches->At(ibranch) ; |
196 | if( fTSTitle.CompareTo(tsMakerBranch->GetTitle())==0 ) |
197 | if( strcmp(tsMakerBranch->GetName(),"AliPHOSTrackSegmentMaker") == 0) |
198 | tsMakerNotFound = kFALSE ; |
199 | } |
200 | if(tsNotFound){ |
201 | tsBranch=(TBranch *) branches->At(ibranch) ; |
202 | if( fTSTitle.CompareTo(tsBranch->GetTitle())==0 ) |
203 | if( strcmp(tsBranch->GetName(),"PHOSTS") == 0) |
204 | tsNotFound = kFALSE ; |
205 | } |
206 | } |
207 | |
208 | if(tsMakerNotFound ||tsNotFound ){ |
209 | cout << "Can't find Branch with TrackSegmentMaker and TrackSegments " ; |
210 | cout << "Do nothing" <<endl ; |
211 | return kFALSE ; |
212 | } |
213 | |
214 | tsMakerBranch->SetAddress(&fTSMaker) ; |
215 | tsBranch->SetAddress(&fTrackSegments) ; |
216 | |
761e34c0 |
217 | tsMakerBranch->GetEntry(0) ; |
218 | tsBranch->GetEntry(0) ; |
7acf6008 |
219 | |
220 | fRecPointsTitle = fTSMaker->GetRecPointsBranch() ; |
1cb7c1ee |
221 | |
7acf6008 |
222 | //reading now recponts branches |
223 | TBranch * emcBranch = 0; |
224 | TBranch * cpvBranch = 0; |
225 | TBranch * cluBranch = 0; |
226 | |
227 | Bool_t emcNotFound = kTRUE ; |
228 | Bool_t cpvNotFound = kTRUE ; |
229 | Bool_t cluNotFound = kTRUE ; |
230 | |
231 | for(ibranch = 0;(ibranch <branches->GetEntries())&&(emcNotFound||cpvNotFound||cluNotFound);ibranch++){ |
232 | if(emcNotFound){ |
233 | emcBranch=(TBranch *) branches->At(ibranch) ; |
234 | if( fRecPointsTitle.CompareTo(emcBranch->GetTitle())==0 ) |
235 | if( strcmp(emcBranch->GetName(),"PHOSEmcRP") == 0) |
236 | emcNotFound = kFALSE ; |
237 | } |
238 | if(cpvNotFound){ |
239 | cpvBranch=(TBranch *) branches->At(ibranch) ; |
240 | if( fRecPointsTitle.CompareTo(cpvBranch->GetTitle())==0 ) |
241 | if( strcmp(cpvBranch->GetName(),"PHOSCpvRP") == 0) |
242 | cpvNotFound = kFALSE ; |
243 | } |
244 | if(cluNotFound){ |
245 | cluBranch=(TBranch *) branches->At(ibranch) ; |
246 | if( fRecPointsTitle.CompareTo(cluBranch->GetTitle())==0 ) |
247 | if( strcmp(cluBranch->GetName(),"AliPHOSClusterizer") == 0) |
248 | cluNotFound = kFALSE ; |
249 | } |
250 | } |
251 | |
252 | if(emcNotFound ||cpvNotFound ||cluNotFound ){ |
253 | cout << "Can't find Branch with RecPoints or AliPHOSClusterizer " ; |
254 | cout << "Do nothing" <<endl ; |
255 | return kFALSE ; |
256 | } |
257 | |
258 | emcBranch->SetAddress(&fEmcRecPoints) ; |
259 | cpvBranch->SetAddress(&fCpvRecPoints) ; |
260 | cluBranch->SetAddress(&fClusterizer) ; |
261 | |
761e34c0 |
262 | emcBranch->GetEntry(0) ; |
263 | cpvBranch->GetEntry(0) ; |
264 | cluBranch->GetEntry(0) ; |
265 | |
7acf6008 |
266 | return kTRUE ; |
267 | |
268 | |
269 | |
270 | } |
69183710 |
271 | //____________________________________________________________________________ |
7acf6008 |
272 | Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const |
69183710 |
273 | { |
274 | // Calculates the distance between the EMC RecPoint and the PPSD RecPoint |
275 | |
69183710 |
276 | TVector3 vecEmc ; |
7acf6008 |
277 | TVector3 vecCpv ; |
278 | |
279 | emc->GetLocalPosition(vecEmc) ; |
280 | cpv->GetLocalPosition(vecCpv) ; |
281 | if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ |
282 | |
283 | // Correct to difference in CPV and EMC position due to different distance to center. |
284 | // we assume, that particle moves from center |
285 | Float_t dCPV = fGeom->GetIPtoOuterCoverDistance(); |
286 | Float_t dEMC = fGeom->GetIPtoCrystalSurface() ; |
287 | dEMC = dEMC / dCPV ; |
288 | vecCpv = dEMC * vecCpv - vecEmc ; |
289 | if (Axis == "X") return vecCpv.X(); |
290 | if (Axis == "Y") return vecCpv.Y(); |
291 | if (Axis == "Z") return vecCpv.Z(); |
292 | if (Axis == "R") return vecCpv.Mag(); |
293 | } |
294 | |
295 | return 100000000 ; |
69183710 |
296 | } |
297 | |
6ad0bfa0 |
298 | //____________________________________________________________________________ |
7acf6008 |
299 | void AliPHOSPIDv1::Exec(Option_t * option) |
6ad0bfa0 |
300 | { |
f035f6ce |
301 | //Steering method |
302 | |
7acf6008 |
303 | if(!fIsInitialized) |
304 | Init() ; |
305 | |
306 | if(strstr(option,"tim")) |
307 | gBenchmark->Start("PHOSPID"); |
308 | |
309 | |
310 | Int_t nEvents = (Int_t) gAlice->TreeE()->GetEntries() ; |
311 | |
312 | for(fNEvent = 0 ;fNEvent <nEvents; fNEvent++){ |
313 | if(!ReadTrackSegments()) |
314 | return ; |
315 | MakeRecParticles() ; |
316 | WriteRecParticles(); |
317 | if(strstr(option,"deb")) |
318 | PrintRecParticles(option) ; |
319 | } |
320 | |
321 | if(strstr(option,"tim")){ |
322 | gBenchmark->Stop("PHOSPID"); |
323 | cout << "AliPHOSPID:" << endl ; |
324 | cout << " took " << gBenchmark->GetCpuTime("PHOSPID") << " seconds for PID " |
325 | << gBenchmark->GetCpuTime("PHOSPID")/nEvents << " seconds per event " << endl ; |
326 | cout << endl ; |
327 | } |
328 | |
329 | } |
330 | //____________________________________________________________________________ |
331 | void AliPHOSPIDv1::MakeRecParticles(){ |
332 | |
b2a60966 |
333 | // Makes a RecParticle out of a TrackSegment |
6ad0bfa0 |
334 | |
7acf6008 |
335 | TIter next(fTrackSegments) ; |
336 | AliPHOSTrackSegment * ts ; |
6ad0bfa0 |
337 | Int_t index = 0 ; |
09fc14a0 |
338 | AliPHOSRecParticle * rp ; |
7acf6008 |
339 | |
340 | Bool_t ellips = fIDOptions.Contains("ell",TString::kIgnoreCase ) ; |
341 | Bool_t disp = fIDOptions.Contains("dis",TString::kIgnoreCase ) ; |
342 | |
343 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { |
fad3e5b9 |
344 | |
7acf6008 |
345 | new( (*fRecParticles)[index] ) AliPHOSRecParticle() ; |
346 | rp = (AliPHOSRecParticle *)fRecParticles->At(index) ; |
347 | rp->SetTraskSegment(index) ; |
f035f6ce |
348 | |
7acf6008 |
349 | AliPHOSEmcRecPoint * emc = 0 ; |
350 | if(ts->GetEmcIndex()>=0) |
351 | emc = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(ts->GetEmcIndex()) ; |
fad3e5b9 |
352 | |
7acf6008 |
353 | AliPHOSRecPoint * cpv = 0 ; |
354 | if(ts->GetCpvIndex()>=0) |
355 | cpv = (AliPHOSRecPoint *) fCpvRecPoints->At(ts->GetCpvIndex()) ; |
fad3e5b9 |
356 | |
7acf6008 |
357 | AliPHOSRecPoint * ppsd = 0 ; |
358 | if(ts->GetPpsdIndex()>=0) |
359 | ppsd= (AliPHOSRecPoint *) fCpvRecPoints->At(ts->GetPpsdIndex()) ; |
360 | |
361 | //set momentum and energy first |
362 | Float_t e = emc->GetEnergy() ; |
363 | TVector3 dir = GetMomentumDirection(emc,cpv,ppsd) ; |
364 | dir.SetMag(e) ; |
365 | |
366 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ; |
367 | rp->SetCalcMass(0); |
368 | |
369 | //now set type (reconstructed) of the particle |
370 | Int_t showerprofile = 0; // 0 narrow and 1 wide |
fad3e5b9 |
371 | |
7acf6008 |
372 | if(ellips){ |
373 | Float_t lambda[2] ; |
374 | emc->GetElipsAxis(lambda) ; |
375 | if(fFormula->Eval(lambda[0],lambda[1]) <= 0 ) |
376 | showerprofile = 1 ; // not narrow |
1cb7c1ee |
377 | } |
fad3e5b9 |
378 | |
7acf6008 |
379 | if(disp) |
380 | if(emc->GetDispersion() > fDispersion ) |
381 | showerprofile = 1 ; // not narrow |
382 | |
383 | |
2aad621e |
384 | // Looking at the photon conversion detector |
7acf6008 |
385 | Int_t pcdetector= 0 ; //1 hit and 0 no hit |
386 | if(ppsd) |
387 | if(GetDistance(emc, ppsd, "R") < fCpvEmcDistance) |
388 | pcdetector = 1 ; |
389 | |
390 | // Looking at the CPV detector |
391 | Int_t cpvdetector= 0 ; //1 hit and 0 no hit |
392 | if(cpv) |
393 | if(GetDistance(emc, cpv, "R") < fCpvEmcDistance) |
394 | cpvdetector = 1 ; |
395 | |
396 | Int_t type = showerprofile + 2 * pcdetector + 4 * cpvdetector ; |
09fc14a0 |
397 | rp->SetType(type) ; |
6ad0bfa0 |
398 | index++ ; |
399 | } |
7acf6008 |
400 | |
6ad0bfa0 |
401 | } |
402 | |
09fc14a0 |
403 | //____________________________________________________________________________ |
7acf6008 |
404 | void AliPHOSPIDv1:: Print(Option_t * option) const |
09fc14a0 |
405 | { |
b2a60966 |
406 | // Print the parameters used for the particle type identification |
f035f6ce |
407 | cout << "=============== AliPHOSPID1 ================" << endl ; |
408 | cout << "Making PID "<< endl ; |
409 | cout << " Headers file: " << fHeaderFileName.Data() << endl ; |
410 | cout << " RecPoints branch title: " << fRecPointsTitle.Data() << endl ; |
411 | cout << " TrackSegments Branch title: " << fTSTitle.Data() << endl ; |
412 | cout << " RecParticles Branch title " << fRecparticlesTitle.Data() << endl; |
413 | cout << "with parameters: " << endl ; |
414 | cout << " Maximal EMC - CPV (PPSD) distance (cm) " << fCpvEmcDistance << endl ; |
415 | if(fIDOptions.Contains("dis",TString::kIgnoreCase )) |
416 | cout << " dispersion cut " << fDispersion << endl ; |
417 | if(fIDOptions.Contains("ell",TString::kIgnoreCase )){ |
418 | cout << "Eliptic cuts function: " << endl ; |
419 | cout << fFormula->GetTitle() << endl ; |
420 | } |
421 | cout << "============================================" << endl ; |
09fc14a0 |
422 | } |
423 | |
424 | //____________________________________________________________________________ |
a4e98857 |
425 | void AliPHOSPIDv1::SetShowerProfileCut(char * formula) |
426 | { |
7acf6008 |
427 | //set shape of the cut on the axis of ellipce, drown around shouer |
428 | //shower considered "narrow" if Formula(lambda[0],lambda[1]) > 0. |
429 | if(fFormula) |
430 | delete fFormula; |
f035f6ce |
431 | fFormula = new TFormula("Lambda Cut",formula) ; |
7acf6008 |
432 | } |
433 | //____________________________________________________________________________ |
434 | void AliPHOSPIDv1::WriteRecParticles() |
09fc14a0 |
435 | { |
7acf6008 |
436 | //check, if these branches already exist |
437 | TBranch * pidBranch = 0; |
438 | TBranch * rpBranch = 0; |
439 | |
440 | TObjArray * branches = gAlice->TreeR()->GetListOfBranches() ; |
441 | Int_t ibranch; |
442 | Bool_t pidNotFound = kTRUE ; |
443 | Bool_t rpNotFound = kTRUE ; |
444 | |
445 | for(ibranch = 0;(ibranch <branches->GetEntries())&& pidNotFound && rpNotFound;ibranch++){ |
446 | if(pidNotFound){ |
447 | pidBranch=(TBranch *) branches->At(ibranch) ; |
448 | if( (strcmp(pidBranch->GetName(),"PHOSPID") == 0) && |
449 | (fRecparticlesTitle.CompareTo(pidBranch->GetTitle()) == 0) ) |
450 | pidNotFound = kFALSE ; |
451 | } |
452 | if(rpNotFound){ |
453 | rpBranch=(TBranch *) branches->At(ibranch) ; |
454 | if( (strcmp(rpBranch->GetName(),"PHOSRP") == 0) && |
455 | (fRecparticlesTitle.CompareTo(rpBranch->GetTitle())==0 )) |
456 | rpNotFound = kFALSE ; |
457 | } |
458 | } |
459 | |
460 | if(!pidNotFound || !rpNotFound) { |
461 | cout << "AliPHOSPIDv1 error: " << endl ; |
462 | cout << " Branch PHOSRP and PHOSPID with title '"<<fRecparticlesTitle.Data()<<"' already exist "<< endl ; |
463 | cout << " can not overwrite " << endl ; |
464 | return ; |
465 | } |
69183710 |
466 | |
7acf6008 |
467 | //Make branch in TreeR for TrackSegments |
468 | char * filename = 0; |
469 | if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){ //generating file name |
470 | filename = new char[strlen(gAlice->GetBaseFile())+20] ; |
471 | sprintf(filename,"%s/PHOS.Reco.root",gAlice->GetBaseFile()) ; |
472 | } |
473 | |
474 | TDirectory *cwd = gDirectory; |
475 | |
476 | //First rp |
477 | Int_t bufferSize = 32000 ; |
478 | rpBranch = gAlice->TreeR()->Branch("PHOSRP",&fRecParticles,bufferSize); |
479 | rpBranch->SetTitle(fRecparticlesTitle.Data()); |
480 | if (filename) { |
481 | rpBranch->SetFile(filename); |
482 | TIter next( rpBranch->GetListOfBranches()); |
761e34c0 |
483 | TBranch * sb ; |
484 | while ((sb=(TBranch*)next())) { |
485 | sb->SetFile(filename); |
7acf6008 |
486 | } |
487 | cwd->cd(); |
488 | } |
489 | |
490 | //second, pid |
491 | Int_t splitlevel = 0 ; |
492 | AliPHOSPIDv1 * pid = this ; |
493 | pidBranch = gAlice->TreeR()->Branch("AliPHOSPID","AliPHOSPIDv1",&pid,bufferSize,splitlevel); |
494 | pidBranch->SetTitle(fRecparticlesTitle.Data()); |
495 | if (filename) { |
496 | pidBranch->SetFile(filename); |
497 | TIter next( pidBranch->GetListOfBranches()); |
761e34c0 |
498 | TBranch * sb ; |
499 | while ((sb=(TBranch*)next())) { |
500 | sb->SetFile(filename); |
7acf6008 |
501 | } |
502 | cwd->cd(); |
503 | } |
504 | |
761e34c0 |
505 | rpBranch->Fill() ; |
506 | pidBranch->Fill() ; |
eec3ac52 |
507 | |
7acf6008 |
508 | gAlice->TreeR()->Write(0,kOverwrite) ; |
509 | |
510 | } |
69183710 |
511 | //____________________________________________________________________________ |
7acf6008 |
512 | void AliPHOSPIDv1::PlotDispersionCuts()const |
69183710 |
513 | { |
a4e98857 |
514 | // produces a plot of the dispersion cut |
515 | TCanvas* lambdas = new TCanvas("lambdas","Cuts on the ellipse axis",200,10,700,500); |
dd5c4038 |
516 | |
7acf6008 |
517 | if(fIDOptions.Contains("ell",TString::kIgnoreCase ) ){ |
518 | TF2 * ell = new TF2("Elliptic Cuts",fFormula->GetName(),0,3,0,3) ; |
519 | ell->SetMinimum(0.0000001) ; |
520 | ell->SetMaximum(0.001) ; |
521 | ell->SetLineStyle(1) ; |
522 | ell->SetLineWidth(2) ; |
523 | ell->Draw() ; |
524 | } |
525 | |
526 | if( fIDOptions.Contains("dis",TString::kIgnoreCase ) ){ |
527 | TF2 * dsp = new TF2("dispersion","(y<x)*(x*x+y*y < [0]*[0])",0,3,0,3) ; |
528 | dsp->SetParameter(0,fDispersion) ; |
529 | dsp->SetMinimum(0.0000001) ; |
530 | dsp->SetMaximum(0.001) ; |
531 | dsp->SetLineStyle(1) ; |
532 | dsp->SetLineColor(2) ; |
533 | dsp->SetLineWidth(2) ; |
534 | dsp->SetNpx(200) ; |
535 | dsp->SetNpy(200) ; |
536 | if(fIDOptions.Contains("ell",TString::kIgnoreCase ) ) |
537 | dsp->Draw("same") ; |
538 | else |
539 | dsp->Draw() ; |
540 | } |
541 | lambdas->Update(); |
542 | } |
69183710 |
543 | |
7acf6008 |
544 | //____________________________________________________________________________ |
545 | TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv,AliPHOSRecPoint * ppsd)const |
546 | { |
547 | // Calculates the momentum direction: |
548 | // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint |
549 | // 2. if a EMC RecPoint and one PPSD RecPoint, direction is given by the line through the 2 recpoints |
550 | // 3. if a EMC RecPoint and two PPSD RecPoints, dirrection is given by the average line through |
551 | // the 2 pairs of recpoints |
552 | // However because of the poor position resolution of PPSD the direction is always taken as if we were |
553 | // in case 1. |
554 | |
555 | TVector3 dir(0,0,0) ; |
556 | |
557 | TVector3 emcglobalpos ; |
558 | TMatrix dummy ; |
559 | |
560 | emc->GetGlobalPosition(emcglobalpos, dummy) ; |
561 | |
69183710 |
562 | |
a4e98857 |
563 | // The following commented code becomes valid once the PPSD provides |
7acf6008 |
564 | // a reasonable position resolution, at least as good as EMC ! |
565 | // TVector3 ppsdlglobalpos ; |
566 | // TVector3 ppsduglobalpos ; |
567 | // if( fPpsdLowRecPoint ){ // certainly a photon that has concerted |
568 | // fPpsdLowRecPoint->GetGlobalPosition(ppsdlglobalpos, mdummy) ; |
569 | // dir = emcglobalpos - ppsdlglobalpos ; |
570 | // if( fPpsdUpRecPoint ){ // not looks like a charged |
571 | // fPpsdUpRecPoint->GetGlobalPosition(ppsduglobalpos, mdummy) ; |
572 | // dir = ( dir + emcglobalpos - ppsduglobalpos ) * 0.5 ; |
573 | // } |
574 | // } |
575 | // else { // looks like a neutral |
576 | // dir = emcglobalpos ; |
577 | // } |
578 | |
579 | dir = emcglobalpos ; |
580 | dir.SetZ( -dir.Z() ) ; // why ? |
581 | dir.SetMag(1.) ; |
582 | |
583 | //account correction to the position of IP |
584 | Float_t xo,yo,zo ; //Coordinates of the origin |
585 | gAlice->Generator()->GetOrigin(xo,yo,zo) ; |
586 | TVector3 origin(xo,yo,zo); |
587 | dir = dir - origin ; |
588 | |
589 | return dir ; |
590 | } |
591 | //____________________________________________________________________________ |
a4e98857 |
592 | void AliPHOSPIDv1::PrintRecParticles(Option_t * option) |
593 | { |
dd5c4038 |
594 | // Print table of reconstructed particles |
595 | |
7acf6008 |
596 | cout << "AliPHOSPIDv1: " << endl ; |
597 | cout << " found " << fRecParticles->GetEntriesFast() << " RecParticles " << endl ; |
598 | |
599 | if(strstr(option,"all")) { // printing found TS |
600 | |
601 | cout << " PARTICLE " |
602 | << " Index " << endl ; |
603 | // << " X " |
604 | // << " Y " |
605 | // << " Z " |
606 | // << " # of primaries " |
607 | // << " Primaries list " << endl; |
608 | |
609 | Int_t index ; |
610 | for (index = 0 ; index < fRecParticles->GetEntries() ; index++) { |
611 | AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ; |
612 | |
613 | Text_t particle[11]; |
614 | switch(rp->GetType()) { |
615 | case AliPHOSFastRecParticle::kNEUTRALEM: |
616 | strcpy( particle, "NEUTRAL_EM"); |
617 | break; |
618 | case AliPHOSFastRecParticle::kNEUTRALHA: |
619 | strcpy(particle, "NEUTRAL_HA"); |
620 | break; |
621 | case AliPHOSFastRecParticle::kGAMMA: |
622 | strcpy(particle, "GAMMA"); |
623 | break ; |
624 | case AliPHOSFastRecParticle::kGAMMAHA: |
625 | strcpy(particle, "GAMMA_H"); |
626 | break ; |
627 | case AliPHOSFastRecParticle::kABSURDEM: |
628 | strcpy(particle, "ABSURD_EM") ; |
629 | break ; |
630 | case AliPHOSFastRecParticle::kABSURDHA: |
631 | strcpy(particle, "ABSURD_HA") ; |
632 | break ; |
633 | case AliPHOSFastRecParticle::kELECTRON: |
634 | strcpy(particle, "ELECTRON") ; |
635 | break ; |
636 | case AliPHOSFastRecParticle::kCHARGEDHA: |
637 | strcpy(particle, "CHARGED_HA") ; |
638 | break ; |
639 | } |
640 | |
641 | // Int_t * primaries; |
642 | // Int_t nprimaries; |
643 | // primaries = rp->GetPrimaries(nprimaries); |
644 | |
645 | cout << setw(15) << particle << " " |
646 | << setw(3) << rp->GetIndexInList() << " " ; |
647 | // << setw(4) << nprimaries << " "; |
648 | // for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) |
649 | // cout << setw(4) << primaries[iprimary] << " "; |
650 | cout << endl; |
651 | } |
652 | cout << "-------------------------------------------" << endl ; |
653 | } |
654 | |
69183710 |
655 | } |
656 | |
7acf6008 |
657 | |
658 | |