Line data Source code
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 : //-----------------------------------------------------//
17 : // //
18 : // Date : March 25 2004 //
19 : // This reads the file PMD.RecPoints.root(TreeR), //
20 : // calls the Clustering algorithm and stores the //
21 : // clustering output in PMD.RecPoints.root(TreeR) //
22 : // //
23 : //-----------------------------------------------------//
24 :
25 : #include <Riostream.h>
26 : #include <TMath.h>
27 : #include <TTree.h>
28 : #include <TObjArray.h>
29 : #include <TClonesArray.h>
30 : #include <TFile.h>
31 : #include <TBranch.h>
32 : #include <TNtuple.h>
33 : #include <TParticle.h>
34 :
35 : #include <TGeoMatrix.h>
36 :
37 : #include "AliGeomManager.h"
38 :
39 : #include "AliPMDcluster.h"
40 : #include "AliPMDclupid.h"
41 : #include "AliPMDrecpoint1.h"
42 : #include "AliPMDrecdata.h"
43 : #include "AliPMDrechit.h"
44 : #include "AliPMDUtility.h"
45 : #include "AliPMDDiscriminator.h"
46 : #include "AliPMDEmpDiscriminator.h"
47 : #include "AliPMDtracker.h"
48 :
49 : #include "AliESDPmdTrack.h"
50 : #include "AliESDEvent.h"
51 : #include "AliLog.h"
52 :
53 12 : ClassImp(AliPMDtracker)
54 :
55 2 : AliPMDtracker::AliPMDtracker():
56 2 : fTreeR(0),
57 6 : fRecpoints(new TClonesArray("AliPMDrecpoint1", 10)),
58 6 : fRechits(new TClonesArray("AliPMDrechit", 10)),
59 6 : fPMDcontin(new TObjArray()),
60 6 : fPMDcontout(new TObjArray()),
61 6 : fPMDutil(new AliPMDUtility()),
62 2 : fPMDrecpoint(0),
63 2 : fPMDclin(0),
64 2 : fPMDclout(0),
65 2 : fXvertex(0.),
66 2 : fYvertex(0.),
67 2 : fZvertex(0.),
68 2 : fSigmaX(0.),
69 2 : fSigmaY(0.),
70 2 : fSigmaZ(0.)
71 10 : {
72 : //
73 : // Default Constructor
74 : //
75 4 : }
76 : //--------------------------------------------------------------------//
77 : AliPMDtracker:: AliPMDtracker(const AliPMDtracker & /* tracker */):
78 0 : TObject(/* tracker */),
79 0 : fTreeR(0),
80 0 : fRecpoints(NULL),
81 0 : fRechits(NULL),
82 0 : fPMDcontin(NULL),
83 0 : fPMDcontout(NULL),
84 0 : fPMDutil(NULL),
85 0 : fPMDrecpoint(0),
86 0 : fPMDclin(0),
87 0 : fPMDclout(0),
88 0 : fXvertex(0.),
89 0 : fYvertex(0.),
90 0 : fZvertex(0.),
91 0 : fSigmaX(0.),
92 0 : fSigmaY(0.),
93 0 : fSigmaZ(0.)
94 0 : {
95 : // copy constructor
96 0 : AliError("Copy constructor not allowed");
97 0 : }
98 :
99 : //--------------------------------------------------------------------//
100 : AliPMDtracker& AliPMDtracker::operator=(const AliPMDtracker & /* tracker */)
101 : {
102 : // assignment operator
103 0 : AliError("Assignment operator not allowed");
104 0 : return *this;
105 : }
106 :
107 : //--------------------------------------------------------------------//
108 : AliPMDtracker::~AliPMDtracker()
109 8 : {
110 : // Destructor
111 2 : if (fRecpoints)
112 : {
113 2 : fRecpoints->Clear();
114 : }
115 2 : if (fRechits)
116 : {
117 2 : fRechits->Clear();
118 : }
119 :
120 2 : if (fPMDcontin)
121 : {
122 2 : fPMDcontin->Delete();
123 4 : delete fPMDcontin;
124 2 : fPMDcontin=0;
125 :
126 2 : }
127 2 : if (fPMDcontout)
128 : {
129 2 : fPMDcontout->Delete();
130 4 : delete fPMDcontout;
131 2 : fPMDcontout=0;
132 :
133 2 : }
134 4 : delete fPMDutil;
135 4 : }
136 : //--------------------------------------------------------------------//
137 : void AliPMDtracker::LoadClusters(TTree *treein)
138 : {
139 : // Load the Reconstructed tree
140 16 : fTreeR = treein;
141 8 : }
142 : //--------------------------------------------------------------------//
143 : void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
144 : {
145 : // Converts digits to recpoints after running clustering
146 : // algorithm on CPV plane and PREshower plane
147 : //
148 :
149 : Int_t idet = 0;
150 : Int_t ismn = 0;
151 16 : Int_t trackno = 1, trackpid = 0;
152 8 : Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
153 :
154 : Int_t *irow;
155 : Int_t *icol;
156 : Int_t *itra;
157 : Int_t *ipid;
158 : Float_t *cadc;
159 :
160 : AliPMDrechit *rechit = 0x0;
161 :
162 8 : TBranch *branch = fTreeR->GetBranch("PMDRecpoint");
163 8 : if (!branch)
164 : {
165 0 : AliError("PMDRecpoint branch not found");
166 0 : return;
167 : }
168 8 : branch->SetAddress(&fRecpoints);
169 :
170 8 : TBranch *branch1 = fTreeR->GetBranch("PMDRechit");
171 8 : if (!branch1)
172 : {
173 0 : AliError("PMDRechit branch not found");
174 0 : return;
175 : }
176 8 : branch1->SetAddress(&fRechits);
177 :
178 : Int_t ncrhit = 0;
179 8 : Int_t nmodules = (Int_t) branch->GetEntries();
180 :
181 24 : AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
182 212 : for (Int_t imodule = 0; imodule < nmodules; imodule++)
183 : {
184 98 : branch->GetEntry(imodule);
185 98 : Int_t nentries = fRecpoints->GetLast();
186 294 : AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
187 : ,nentries));
188 :
189 540 : for(Int_t ient = 0; ient < nentries+1; ient++)
190 : {
191 172 : fPMDrecpoint = (AliPMDrecpoint1*)fRecpoints->UncheckedAt(ient);
192 172 : idet = fPMDrecpoint->GetDetector();
193 172 : ismn = fPMDrecpoint->GetSMNumber();
194 172 : clusdata[0] = fPMDrecpoint->GetClusX();
195 172 : clusdata[1] = fPMDrecpoint->GetClusY();
196 172 : clusdata[2] = fPMDrecpoint->GetClusADC();
197 172 : clusdata[3] = fPMDrecpoint->GetClusCells();
198 172 : clusdata[4] = fPMDrecpoint->GetClusSigmaX();
199 172 : clusdata[5] = fPMDrecpoint->GetClusSigmaY();
200 :
201 344 : if (clusdata[4] >= 0. && clusdata[5] >= 0.)
202 : {
203 : // extract the associated cell information
204 172 : branch1->GetEntry(ncrhit);
205 172 : Int_t nenbr1 = fRechits->GetLast() + 1;
206 :
207 172 : irow = new Int_t[nenbr1];
208 172 : icol = new Int_t[nenbr1];
209 172 : itra = new Int_t[nenbr1];
210 172 : ipid = new Int_t[nenbr1];
211 172 : cadc = new Float_t[nenbr1];
212 :
213 1048 : for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
214 : {
215 352 : irow[ient1] = -99;
216 352 : icol[ient1] = -99;
217 352 : itra[ient1] = -99;
218 352 : ipid[ient1] = -99;
219 352 : cadc[ient1] = 0.;
220 : }
221 1048 : for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
222 : {
223 352 : rechit = (AliPMDrechit*)fRechits->UncheckedAt(ient1);
224 : //irow[ient1] = rechit->GetCellX();
225 : //icol[ient1] = rechit->GetCellY();
226 352 : itra[ient1] = rechit->GetCellTrack();
227 352 : ipid[ient1] = rechit->GetCellPid();
228 352 : cadc[ient1] = rechit->GetCellAdc();
229 : }
230 172 : if (idet == 0)
231 : {
232 88 : AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
233 : trackno, trackpid);
234 88 : }
235 84 : else if (idet == 1)
236 : {
237 84 : trackno = itra[0];
238 84 : trackpid = ipid[0];
239 84 : }
240 :
241 344 : delete [] irow;
242 344 : delete [] icol;
243 344 : delete [] itra;
244 344 : delete [] ipid;
245 344 : delete [] cadc;
246 :
247 344 : fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
248 172 : fPMDcontin->Add(fPMDclin);
249 :
250 172 : ncrhit++;
251 172 : }
252 : }
253 : }
254 :
255 8 : AliPMDEmpDiscriminator pmddiscriminator;
256 8 : pmddiscriminator.Discrimination(fPMDcontin,fPMDcontout);
257 :
258 : // alignment implemention
259 :
260 8 : Double_t sectr[4][3] = { {0.,0.,0.},{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
261 8 : TString snsector="PMD/Sector";
262 8 : TString symname;
263 8 : TGeoHMatrix gpmdor;
264 :
265 80 : for(Int_t isector=1; isector<=4; isector++)
266 : {
267 32 : symname = snsector;
268 32 : symname += isector;
269 64 : TGeoHMatrix *gpmdal = AliGeomManager::GetMatrix(symname);
270 32 : Double_t *tral = gpmdal->GetTranslation();
271 :
272 64 : AliGeomManager::GetOrigGlobalMatrix(symname, gpmdor);
273 32 : Double_t *tror = gpmdor.GetTranslation();
274 :
275 256 : for(Int_t ixyz=0; ixyz<3; ixyz++)
276 : {
277 96 : sectr[isector-1][ixyz] = tral[ixyz] - tror[ixyz];
278 : }
279 : }
280 :
281 : const Float_t kzpos = 361.5; // middle of the PMD
282 :
283 : Int_t ix = -1, iy = -1;
284 : Int_t det = 0, smn = 0, trno = 1, trpid = 0, mstat = 0;
285 : Float_t xpos = 0., ypos = 0.;
286 : Float_t adc = 0., ncell = 0., radx = 0., rady = 0.;
287 8 : Float_t xglobal = 0., yglobal = 0., zglobal = 0;
288 : Float_t pid = 0.;
289 :
290 8 : fPMDutil->ApplyAlignment(sectr);
291 :
292 8 : Int_t nentries2 = fPMDcontout->GetEntries();
293 40 : AliDebug(1,Form("Number of clusters coming after discrimination = %d"
294 : ,nentries2));
295 360 : for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
296 : {
297 172 : fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
298 :
299 172 : det = fPMDclout->GetDetector();
300 172 : smn = fPMDclout->GetSMN();
301 172 : trno = fPMDclout->GetClusTrackNo();
302 172 : trpid = fPMDclout->GetClusTrackPid();
303 172 : mstat = fPMDclout->GetClusMatching();
304 172 : xpos = fPMDclout->GetClusX();
305 172 : ypos = fPMDclout->GetClusY();
306 172 : adc = fPMDclout->GetClusADC();
307 172 : ncell = fPMDclout->GetClusCells();
308 172 : radx = fPMDclout->GetClusSigmaX();
309 : // Here in the variable "rady" we are keeping the row and col
310 : // of the single isolated cluster having ncell = 1 for offline
311 : // calibration
312 :
313 294 : if ((radx > 999. && radx < 1000.) && ncell == 1)
314 : {
315 122 : if (smn < 12)
316 : {
317 90 : ix = (Int_t) (ypos +0.5);
318 90 : iy = (Int_t) xpos;
319 90 : }
320 32 : else if (smn > 12 && smn < 24)
321 : {
322 28 : ix = (Int_t) xpos;
323 28 : iy = (Int_t) (ypos +0.5);
324 28 : }
325 122 : rady = (Float_t) (ix*100 + iy);
326 122 : }
327 : else
328 : {
329 50 : rady = fPMDclout->GetClusSigmaY();
330 : }
331 172 : pid = fPMDclout->GetClusPID();
332 :
333 : //
334 : /**********************************************************************
335 : * det : Detector, 0: PRE & 1:CPV *
336 : * smn : Serial Module Number 0 to 23 for each plane *
337 : * xpos : x-position of the cluster *
338 : * ypos : y-position of the cluster *
339 : * THESE xpos & ypos are not the true xpos and ypos *
340 : * for some of the unit modules. They are rotated. *
341 : * adc : ADC contained in the cluster *
342 : * ncell : Number of cells contained in the cluster *
343 : * rad : radius of the cluster (1d fit) *
344 : **********************************************************************/
345 : //
346 :
347 :
348 172 : if (det == 0)
349 : {
350 88 : zglobal = kzpos + 1.65; // PREshower plane
351 88 : }
352 84 : else if (det == 1)
353 : {
354 84 : zglobal = kzpos - 1.65; // CPV plane
355 84 : }
356 :
357 172 : fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal,zglobal);
358 :
359 : // Fill ESD
360 :
361 344 : AliESDPmdTrack *esdpmdtr = new AliESDPmdTrack();
362 :
363 172 : esdpmdtr->SetDetector(det);
364 172 : esdpmdtr->SetSmn(smn);
365 172 : esdpmdtr->SetClusterTrackNo(trno);
366 172 : esdpmdtr->SetClusterTrackPid(trpid);
367 172 : esdpmdtr->SetClusterMatching(mstat);
368 :
369 172 : esdpmdtr->SetClusterX(xglobal);
370 172 : esdpmdtr->SetClusterY(yglobal);
371 172 : esdpmdtr->SetClusterZ(zglobal);
372 172 : esdpmdtr->SetClusterADC(adc);
373 172 : esdpmdtr->SetClusterCells(ncell);
374 172 : esdpmdtr->SetClusterPID(pid);
375 172 : esdpmdtr->SetClusterSigmaX(radx);
376 172 : esdpmdtr->SetClusterSigmaY(rady);
377 :
378 172 : event->AddPmdTrack(esdpmdtr);
379 344 : delete esdpmdtr;
380 : }
381 :
382 8 : fPMDcontin->Delete();
383 8 : fPMDcontout->Delete();
384 :
385 16 : }
386 : //--------------------------------------------------------------------//
387 : void AliPMDtracker::AssignTrPidToCluster(Int_t nentry, Int_t *itra,
388 : Int_t *ipid, Float_t *cadc,
389 : Int_t &trackno, Int_t &trackpid)
390 : {
391 : // assign the track number and the corresponding pid to a cluster
392 : // split cluster part will be done at the time of calculating eff/pur
393 :
394 176 : Int_t *phentry = new Int_t [nentry];
395 88 : Int_t *hadentry = new Int_t [nentry];
396 : Int_t *trpid = 0x0;
397 : Int_t *sortcoord = 0x0;
398 : Float_t *trenergy = 0x0;
399 :
400 : Int_t ngtrack = 0;
401 : Int_t nhtrack = 0;
402 644 : for (Int_t i = 0; i < nentry; i++)
403 : {
404 234 : phentry[i] = -1;
405 234 : hadentry[i] = -1;
406 :
407 234 : if (ipid[i] == 22)
408 : {
409 60 : phentry[ngtrack] = i;
410 60 : ngtrack++;
411 60 : }
412 174 : else if (ipid[i] != 22)
413 : {
414 174 : hadentry[nhtrack] = i;
415 174 : nhtrack++;
416 174 : }
417 : }
418 :
419 88 : Int_t nghadtrack = ngtrack + nhtrack;
420 :
421 88 : if (ngtrack == 0)
422 : {
423 : // hadron track
424 : // no need of track number, set to -1
425 75 : trackpid = 8;
426 75 : trackno = -1;
427 75 : }
428 13 : else if (ngtrack >= 1)
429 : {
430 : // one or more than one photon track + charged track
431 : // find out which track deposits maximum energy and
432 : // assign that track number and track pid
433 :
434 13 : trenergy = new Float_t [nghadtrack];
435 13 : trpid = new Int_t [nghadtrack];
436 13 : sortcoord = new Int_t [2*nghadtrack];
437 146 : for (Int_t i = 0; i < ngtrack; i++)
438 : {
439 60 : trenergy[i] = 0.;
440 60 : trpid[i] = -1;
441 1524 : for (Int_t j = 0; j < nentry; j++)
442 : {
443 1404 : if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
444 : {
445 702 : trenergy[i] += cadc[j];
446 702 : trpid[i] = 22;
447 702 : }
448 : }
449 : }
450 26 : for (Int_t i = ngtrack; i < nghadtrack; i++)
451 : {
452 0 : trenergy[i] = 0.;
453 0 : trpid[i] = -1;
454 0 : for (Int_t j = 0; j < nentry; j++)
455 : {
456 0 : if (ipid[j] != 22 && itra[j] == itra[hadentry[i-ngtrack]])
457 : {
458 0 : trenergy[i] += cadc[j];
459 0 : trpid[i] = ipid[j];
460 0 : }
461 : }
462 : }
463 :
464 : Bool_t jsort = true;
465 13 : TMath::Sort(nghadtrack,trenergy,sortcoord,jsort);
466 :
467 13 : Int_t gtr = sortcoord[0];
468 13 : if (trpid[gtr] == 22)
469 : {
470 13 : trackpid = 22;
471 13 : trackno = itra[phentry[gtr]]; // highest adc track
472 13 : }
473 : else
474 : {
475 0 : trackpid = 8;
476 0 : trackno = -1;
477 : }
478 :
479 26 : delete [] trenergy;
480 26 : delete [] trpid;
481 26 : delete [] sortcoord;
482 :
483 13 : } // end of ngtrack >= 1
484 :
485 176 : delete [] phentry;
486 176 : delete [] hadentry;
487 :
488 88 : }
489 : //--------------------------------------------------------------------//
490 : void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
491 : {
492 0 : fXvertex = vtx[0];
493 0 : fYvertex = vtx[1];
494 0 : fZvertex = vtx[2];
495 0 : fSigmaX = evtx[0];
496 0 : fSigmaY = evtx[1];
497 0 : fSigmaZ = evtx[2];
498 0 : }
499 : //--------------------------------------------------------------------//
500 : void AliPMDtracker::ResetClusters()
501 : {
502 0 : if (fRecpoints) fRecpoints->Clear();
503 0 : }
504 : //--------------------------------------------------------------------//
|
