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