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 <TGeoMatrix.h>
37 #include "AliGeomManager.h"
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"
49 #include "AliESDPmdTrack.h"
50 #include "AliESDEvent.h"
53 ClassImp(AliPMDtracker)
55 AliPMDtracker::AliPMDtracker():
57 fRecpoints(new TClonesArray("AliPMDrecpoint1", 10)),
58 fRechits(new TClonesArray("AliPMDrechit", 10)),
59 fPMDcontin(new TObjArray()),
60 fPMDcontout(new TObjArray()),
61 fPMDutil(new AliPMDUtility()),
73 // Default Constructor
76 //--------------------------------------------------------------------//
77 AliPMDtracker:: AliPMDtracker(const AliPMDtracker & /* tracker */):
78 TObject(/* tracker */),
96 AliError("Copy constructor not allowed");
99 //--------------------------------------------------------------------//
100 AliPMDtracker& AliPMDtracker::operator=(const AliPMDtracker & /* tracker */)
102 // assignment operator
103 AliError("Assignment operator not allowed");
107 //--------------------------------------------------------------------//
108 AliPMDtracker::~AliPMDtracker()
122 fPMDcontin->Delete();
129 fPMDcontout->Delete();
136 //--------------------------------------------------------------------//
137 void AliPMDtracker::LoadClusters(TTree *treein)
139 // Load the Reconstructed tree
142 //--------------------------------------------------------------------//
143 void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
145 // Converts digits to recpoints after running clustering
146 // algorithm on CPV plane and PREshower plane
151 Int_t trackno = 1, trackpid = 0;
152 Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
160 AliPMDrechit *rechit = 0x0;
162 TBranch *branch = fTreeR->GetBranch("PMDRecpoint");
165 AliError("PMDRecpoint branch not found");
168 branch->SetAddress(&fRecpoints);
170 TBranch *branch1 = fTreeR->GetBranch("PMDRechit");
173 AliError("PMDRechit branch not found");
176 branch1->SetAddress(&fRechits);
179 Int_t nmodules = (Int_t) branch->GetEntries();
181 AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
182 for (Int_t imodule = 0; imodule < nmodules; imodule++)
184 branch->GetEntry(imodule);
185 Int_t nentries = fRecpoints->GetLast();
186 AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
189 for(Int_t ient = 0; ient < nentries+1; ient++)
191 fPMDrecpoint = (AliPMDrecpoint1*)fRecpoints->UncheckedAt(ient);
192 idet = fPMDrecpoint->GetDetector();
193 ismn = fPMDrecpoint->GetSMNumber();
194 clusdata[0] = fPMDrecpoint->GetClusX();
195 clusdata[1] = fPMDrecpoint->GetClusY();
196 clusdata[2] = fPMDrecpoint->GetClusADC();
197 clusdata[3] = fPMDrecpoint->GetClusCells();
198 clusdata[4] = fPMDrecpoint->GetClusSigmaX();
199 clusdata[5] = fPMDrecpoint->GetClusSigmaY();
201 if (clusdata[4] >= 0. && clusdata[5] >= 0.)
203 // extract the associated cell information
204 branch1->GetEntry(ncrhit);
205 Int_t nenbr1 = fRechits->GetLast() + 1;
207 irow = new Int_t[nenbr1];
208 icol = new Int_t[nenbr1];
209 itra = new Int_t[nenbr1];
210 ipid = new Int_t[nenbr1];
211 cadc = new Float_t[nenbr1];
213 for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
221 for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
223 rechit = (AliPMDrechit*)fRechits->UncheckedAt(ient1);
224 //irow[ient1] = rechit->GetCellX();
225 //icol[ient1] = rechit->GetCellY();
226 itra[ient1] = rechit->GetCellTrack();
227 ipid[ient1] = rechit->GetCellPid();
228 cadc[ient1] = rechit->GetCellAdc();
232 AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
247 fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
248 fPMDcontin->Add(fPMDclin);
255 AliPMDEmpDiscriminator pmddiscriminator;
256 pmddiscriminator.Discrimination(fPMDcontin,fPMDcontout);
258 // alignment implemention
260 Double_t sectr[4][3] = { {0.,0.,0.},{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
261 TString snsector="PMD/Sector";
265 for(Int_t isector=1; isector<=4; isector++)
269 TGeoHMatrix *gpmdal = AliGeomManager::GetMatrix(symname);
270 Double_t *tral = gpmdal->GetTranslation();
272 AliGeomManager::GetOrigGlobalMatrix(symname, gpmdor);
273 Double_t *tror = gpmdor.GetTranslation();
275 for(Int_t ixyz=0; ixyz<3; ixyz++)
277 sectr[isector-1][ixyz] = tral[ixyz] - tror[ixyz];
281 const Float_t kzpos = 361.5; // middle of the PMD
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 Float_t xglobal = 0., yglobal = 0., zglobal = 0;
290 fPMDutil->ApplyAlignment(sectr);
292 Int_t nentries2 = fPMDcontout->GetEntries();
293 AliDebug(1,Form("Number of clusters coming after discrimination = %d"
295 for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
297 fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
299 det = fPMDclout->GetDetector();
300 smn = fPMDclout->GetSMN();
301 trno = fPMDclout->GetClusTrackNo();
302 trpid = fPMDclout->GetClusTrackPid();
303 mstat = fPMDclout->GetClusMatching();
304 xpos = fPMDclout->GetClusX();
305 ypos = fPMDclout->GetClusY();
306 adc = fPMDclout->GetClusADC();
307 ncell = fPMDclout->GetClusCells();
308 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
313 if ((radx > 999. && radx < 1000.) && ncell == 1)
317 ix = (Int_t) (ypos +0.5);
320 else if (smn > 12 && smn < 24)
323 iy = (Int_t) (ypos +0.5);
325 rady = (Float_t) (ix*100 + iy);
329 rady = fPMDclout->GetClusSigmaY();
331 pid = fPMDclout->GetClusPID();
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 **********************************************************************/
350 zglobal = kzpos + 1.65; // PREshower plane
354 zglobal = kzpos - 1.65; // CPV plane
357 fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal,zglobal);
361 AliESDPmdTrack *esdpmdtr = new AliESDPmdTrack();
363 esdpmdtr->SetDetector(det);
364 esdpmdtr->SetSmn(smn);
365 esdpmdtr->SetClusterTrackNo(trno);
366 esdpmdtr->SetClusterTrackPid(trpid);
367 esdpmdtr->SetClusterMatching(mstat);
369 esdpmdtr->SetClusterX(xglobal);
370 esdpmdtr->SetClusterY(yglobal);
371 esdpmdtr->SetClusterZ(zglobal);
372 esdpmdtr->SetClusterADC(adc);
373 esdpmdtr->SetClusterCells(ncell);
374 esdpmdtr->SetClusterPID(pid);
375 esdpmdtr->SetClusterSigmaX(radx);
376 esdpmdtr->SetClusterSigmaY(rady);
378 event->AddPmdTrack(esdpmdtr);
382 fPMDcontin->Delete();
383 fPMDcontout->Delete();
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)
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
394 Int_t *phentry = new Int_t [nentry];
395 Int_t *hadentry = new Int_t [nentry];
397 Int_t *sortcoord = 0x0;
398 Float_t *trenergy = 0x0;
402 for (Int_t i = 0; i < nentry; i++)
409 phentry[ngtrack] = i;
412 else if (ipid[i] != 22)
414 hadentry[nhtrack] = i;
419 Int_t nghadtrack = ngtrack + nhtrack;
424 // no need of track number, set to -1
428 else if (ngtrack >= 1)
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
434 trenergy = new Float_t [nghadtrack];
435 trpid = new Int_t [nghadtrack];
436 sortcoord = new Int_t [2*nghadtrack];
437 for (Int_t i = 0; i < ngtrack; i++)
441 for (Int_t j = 0; j < nentry; j++)
443 if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
445 trenergy[i] += cadc[j];
450 for (Int_t i = ngtrack; i < nghadtrack; i++)
454 for (Int_t j = 0; j < nentry; j++)
456 if (ipid[j] != 22 && itra[j] == itra[hadentry[i-ngtrack]])
458 trenergy[i] += cadc[j];
465 TMath::Sort(nghadtrack,trenergy,sortcoord,jsort);
467 Int_t gtr = sortcoord[0];
468 if (trpid[gtr] == 22)
471 trackno = itra[phentry[gtr]]; // highest adc track
483 } // end of ngtrack >= 1
489 //--------------------------------------------------------------------//
490 void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
499 //--------------------------------------------------------------------//
500 void AliPMDtracker::ResetClusters()
502 if (fRecpoints) fRecpoints->Clear();
504 //--------------------------------------------------------------------//