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 //-----------------------------------------------------//
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) //
23 //-----------------------------------------------------//
25 #include <Riostream.h>
28 #include <TObjArray.h>
29 #include <TClonesArray.h>
33 #include <TParticle.h>
35 #include "AliPMDcluster.h"
36 #include "AliPMDclupid.h"
37 #include "AliPMDrecpoint1.h"
38 #include "AliPMDrecdata.h"
39 #include "AliPMDrechit.h"
40 #include "AliPMDUtility.h"
41 #include "AliPMDDiscriminator.h"
42 #include "AliPMDEmpDiscriminator.h"
43 #include "AliPMDtracker.h"
45 #include "AliESDPmdTrack.h"
46 #include "AliESDEvent.h"
49 ClassImp(AliPMDtracker)
51 AliPMDtracker::AliPMDtracker():
53 fRecpoints(new TClonesArray("AliPMDrecpoint1", 10)),
54 fRechits(new TClonesArray("AliPMDrechit", 10)),
55 fPMDcontin(new TObjArray()),
56 fPMDcontout(new TObjArray()),
57 fPMDutil(new AliPMDUtility()),
69 // Default Constructor
72 //--------------------------------------------------------------------//
73 AliPMDtracker:: AliPMDtracker(const AliPMDtracker & /* tracker */):
74 TObject(/* tracker */),
92 AliError("Copy constructor not allowed");
95 //--------------------------------------------------------------------//
96 AliPMDtracker& AliPMDtracker::operator=(const AliPMDtracker & /* tracker */)
98 // assignment operator
99 AliError("Assignment operator not allowed");
103 //--------------------------------------------------------------------//
104 AliPMDtracker::~AliPMDtracker()
118 fPMDcontin->Delete();
125 fPMDcontout->Delete();
132 //--------------------------------------------------------------------//
133 void AliPMDtracker::LoadClusters(TTree *treein)
135 // Load the Reconstructed tree
138 //--------------------------------------------------------------------//
139 void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
141 // Converts digits to recpoints after running clustering
142 // algorithm on CPV plane and PREshower plane
147 Int_t trackno, trackpid;
156 AliPMDrechit *rechit = 0x0;
158 TBranch *branch = fTreeR->GetBranch("PMDRecpoint");
161 AliError("PMDRecpoint branch not found");
164 branch->SetAddress(&fRecpoints);
166 TBranch *branch1 = fTreeR->GetBranch("PMDRechit");
169 AliError("PMDRechit branch not found");
172 branch1->SetAddress(&fRechits);
175 Int_t nmodules = (Int_t) branch->GetEntries();
177 AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
178 for (Int_t imodule = 0; imodule < nmodules; imodule++)
180 branch->GetEntry(imodule);
181 Int_t nentries = fRecpoints->GetLast();
182 AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
185 for(Int_t ient = 0; ient < nentries+1; ient++)
187 fPMDrecpoint = (AliPMDrecpoint1*)fRecpoints->UncheckedAt(ient);
188 idet = fPMDrecpoint->GetDetector();
189 ismn = fPMDrecpoint->GetSMNumber();
190 clusdata[0] = fPMDrecpoint->GetClusX();
191 clusdata[1] = fPMDrecpoint->GetClusY();
192 clusdata[2] = fPMDrecpoint->GetClusADC();
193 clusdata[3] = fPMDrecpoint->GetClusCells();
194 clusdata[4] = fPMDrecpoint->GetClusSigmaX();
195 clusdata[5] = fPMDrecpoint->GetClusSigmaY();
197 if (clusdata[4] >= 0. && clusdata[5] >= 0.)
199 // extract the associated cell information
200 branch1->GetEntry(ncrhit);
201 Int_t nenbr1 = fRechits->GetLast() + 1;
203 irow = new Int_t[nenbr1];
204 icol = new Int_t[nenbr1];
205 itra = new Int_t[nenbr1];
206 ipid = new Int_t[nenbr1];
207 cadc = new Float_t[nenbr1];
209 for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
211 rechit = (AliPMDrechit*)fRechits->UncheckedAt(ient1);
212 //irow[ient1] = rechit->GetCellX();
213 //icol[ient1] = rechit->GetCellY();
214 itra[ient1] = rechit->GetCellTrack();
215 ipid[ient1] = rechit->GetCellPid();
216 cadc[ient1] = rechit->GetCellAdc();
220 AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
235 fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
236 fPMDcontin->Add(fPMDclin);
243 AliPMDEmpDiscriminator pmddiscriminator;
244 pmddiscriminator.Discrimination(fPMDcontin,fPMDcontout);
246 const Float_t kzpos = 361.5; // middle of the PMD
248 Int_t ix = -1, iy = -1;
249 Int_t det,smn,trno,trpid,mstat;
251 Float_t adc, ncell, radx, rady;
252 Float_t xglobal = 0., yglobal = 0., zglobal = 0;
255 fPMDutil->ApplyAlignment();
257 Int_t nentries2 = fPMDcontout->GetEntries();
258 AliDebug(1,Form("Number of clusters coming after discrimination = %d"
260 for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
262 fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
264 det = fPMDclout->GetDetector();
265 smn = fPMDclout->GetSMN();
266 trno = fPMDclout->GetClusTrackNo();
267 trpid = fPMDclout->GetClusTrackPid();
268 mstat = fPMDclout->GetClusMatching();
269 xpos = fPMDclout->GetClusX();
270 ypos = fPMDclout->GetClusY();
271 adc = fPMDclout->GetClusADC();
272 ncell = fPMDclout->GetClusCells();
273 radx = fPMDclout->GetClusSigmaX();
274 // Here in the variable "rady" we are keeping the row and col
275 // of the single isolated cluster having ncell = 1 for offline
278 if ((radx > 999. && radx < 1000.) && ncell == 1)
282 ix = (Int_t) (ypos +0.5);
285 else if (smn > 12 && smn < 24)
288 iy = (Int_t) (ypos +0.5);
290 rady = (Float_t) (ix*100 + iy);
294 rady = fPMDclout->GetClusSigmaY();
296 pid = fPMDclout->GetClusPID();
299 /**********************************************************************
300 * det : Detector, 0: PRE & 1:CPV *
301 * smn : Serial Module Number 0 to 23 for each plane *
302 * xpos : x-position of the cluster *
303 * ypos : y-position of the cluster *
304 * THESE xpos & ypos are not the true xpos and ypos *
305 * for some of the unit modules. They are rotated. *
306 * adc : ADC contained in the cluster *
307 * ncell : Number of cells contained in the cluster *
308 * rad : radius of the cluster (1d fit) *
309 **********************************************************************/
315 zglobal = kzpos + 1.65; // PREshower plane
319 zglobal = kzpos - 1.65; // CPV plane
322 fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal,zglobal);
326 AliESDPmdTrack *esdpmdtr = new AliESDPmdTrack();
328 esdpmdtr->SetDetector(det);
329 esdpmdtr->SetSmn(smn);
330 esdpmdtr->SetClusterTrackNo(trno);
331 esdpmdtr->SetClusterTrackPid(trpid);
332 esdpmdtr->SetClusterMatching(mstat);
334 esdpmdtr->SetClusterX(xglobal);
335 esdpmdtr->SetClusterY(yglobal);
336 esdpmdtr->SetClusterZ(zglobal);
337 esdpmdtr->SetClusterADC(adc);
338 esdpmdtr->SetClusterCells(ncell);
339 esdpmdtr->SetClusterPID(pid);
340 esdpmdtr->SetClusterSigmaX(radx);
341 esdpmdtr->SetClusterSigmaY(rady);
343 event->AddPmdTrack(esdpmdtr);
347 fPMDcontin->Delete();
348 fPMDcontout->Delete();
351 //--------------------------------------------------------------------//
352 void AliPMDtracker::AssignTrPidToCluster(Int_t nentry, Int_t *itra,
353 Int_t *ipid, Float_t *cadc,
354 Int_t &trackno, Int_t &trackpid)
356 // assign the track number and the corresponding pid to a cluster
357 // split cluster part will be done at the time of calculating eff/pur
359 Int_t *phentry = new Int_t [nentry];
360 Int_t *hadentry = new Int_t [nentry];
361 Int_t *trenergy = 0x0;
363 Int_t *sortcoord = 0x0;
367 for (Int_t i = 0; i < nentry; i++)
374 phentry[ngtrack] = i;
377 else if (ipid[i] != 22)
379 hadentry[nhtrack] = i;
384 Int_t nghadtrack = ngtrack + nhtrack;
389 // no need of track number, set to -1
393 else if (ngtrack >= 1)
395 // one or more than one photon track + charged track
396 // find out which track deposits maximum energy and
397 // assign that track number and track pid
399 trenergy = new Int_t [nghadtrack];
400 trpid = new Int_t [nghadtrack];
401 sortcoord = new Int_t [nghadtrack];
402 for (Int_t i = 0; i < ngtrack; i++)
406 for (Int_t j = 0; j < nentry; j++)
408 if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
410 trenergy[i] += cadc[j];
415 for (Int_t i = ngtrack; i < nghadtrack; i++)
419 for (Int_t j = 0; j < nentry; j++)
421 if (ipid[j] != 22 && itra[j] == itra[hadentry[i-ngtrack]])
423 trenergy[i] += cadc[j];
430 TMath::Sort(nghadtrack,trenergy,sortcoord,jsort);
432 Int_t gtr = sortcoord[0];
433 if (trpid[gtr] == 22)
436 trackno = itra[phentry[gtr]]; // highest adc track
448 } // end of ngtrack >= 1
451 //--------------------------------------------------------------------//
452 void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
461 //--------------------------------------------------------------------//
462 void AliPMDtracker::ResetClusters()
464 if (fRecpoints) fRecpoints->Clear();
466 //--------------------------------------------------------------------//