1. New RCU format
[u/mrichter/AliRoot.git] / PMD / AliPMDtracker.cxx
CommitLineData
7dcaf913 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>
7dcaf913 27#include <TTree.h>
7dcaf913 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 "AliPMDcluster.h"
36#include "AliPMDclupid.h"
37#include "AliPMDrecpoint1.h"
b0e4d1e1 38#include "AliPMDrecdata.h"
39#include "AliPMDrechit.h"
7dcaf913 40#include "AliPMDUtility.h"
41#include "AliPMDDiscriminator.h"
26f1ae17 42#include "AliPMDEmpDiscriminator.h"
7dcaf913 43#include "AliPMDtracker.h"
44
45#include "AliESDPmdTrack.h"
af885e0f 46#include "AliESDEvent.h"
ecee2a1a 47#include "AliLog.h"
7dcaf913 48
49ClassImp(AliPMDtracker)
50
51AliPMDtracker::AliPMDtracker():
52 fTreeR(0),
722ccc67 53 fRecpoints(new TClonesArray("AliPMDrecpoint1", 10)),
b0e4d1e1 54 fRechits(new TClonesArray("AliPMDrechit", 10)),
7dcaf913 55 fPMDcontin(new TObjArray()),
56 fPMDcontout(new TObjArray()),
7dcaf913 57 fPMDutil(new AliPMDUtility()),
58 fPMDrecpoint(0),
59 fPMDclin(0),
60 fPMDclout(0),
7dcaf913 61 fXvertex(0.),
62 fYvertex(0.),
63 fZvertex(0.),
64 fSigmaX(0.),
65 fSigmaY(0.),
66 fSigmaZ(0.)
67{
68 //
69 // Default Constructor
70 //
71}
72//--------------------------------------------------------------------//
a48edddd 73AliPMDtracker:: AliPMDtracker(const AliPMDtracker & /* tracker */):
2332574a 74 TObject(/* tracker */),
75 fTreeR(0),
76 fRecpoints(NULL),
b0e4d1e1 77 fRechits(NULL),
2332574a 78 fPMDcontin(NULL),
79 fPMDcontout(NULL),
80 fPMDutil(NULL),
81 fPMDrecpoint(0),
82 fPMDclin(0),
83 fPMDclout(0),
84 fXvertex(0.),
85 fYvertex(0.),
86 fZvertex(0.),
87 fSigmaX(0.),
88 fSigmaY(0.),
89 fSigmaZ(0.)
a48edddd 90{
91 // copy constructor
92 AliError("Copy constructor not allowed");
93}
94
95//--------------------------------------------------------------------//
96AliPMDtracker& AliPMDtracker::operator=(const AliPMDtracker & /* tracker */)
97{
98 // assignment operator
99 AliError("Assignment operator not allowed");
100 return *this;
101}
102
103//--------------------------------------------------------------------//
7dcaf913 104AliPMDtracker::~AliPMDtracker()
105{
106 // Destructor
107 if (fRecpoints)
108 {
722ccc67 109 fRecpoints->Clear();
7dcaf913 110 }
b0e4d1e1 111 if (fRechits)
112 {
113 fRechits->Clear();
114 }
115
7dcaf913 116 if (fPMDcontin)
117 {
118 fPMDcontin->Delete();
119 delete fPMDcontin;
120 fPMDcontin=0;
a1e9ac09 121
7dcaf913 122 }
123 if (fPMDcontout)
722ccc67 124 {
7dcaf913 125 fPMDcontout->Delete();
126 delete fPMDcontout;
127 fPMDcontout=0;
a1e9ac09 128
7dcaf913 129 }
316c6cd9 130 delete fPMDutil;
7dcaf913 131}
132//--------------------------------------------------------------------//
133void AliPMDtracker::LoadClusters(TTree *treein)
134{
135 // Load the Reconstructed tree
136 fTreeR = treein;
137}
138//--------------------------------------------------------------------//
af885e0f 139void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
7dcaf913 140{
141 // Converts digits to recpoints after running clustering
142 // algorithm on CPV plane and PREshower plane
143 //
144
145 Int_t idet;
146 Int_t ismn;
b0e4d1e1 147 Int_t trackno, trackpid;
26f1ae17 148 Float_t clusdata[6];
b0e4d1e1 149
150 Int_t *irow;
151 Int_t *icol;
152 Int_t *itra;
153 Int_t *ipid;
154 Float_t *cadc;
155
156 AliPMDrechit *rechit = 0x0;
7dcaf913 157
158 TBranch *branch = fTreeR->GetBranch("PMDRecpoint");
ecee2a1a 159 if (!branch)
160 {
161 AliError("PMDRecpoint branch not found");
162 return;
163 }
7dcaf913 164 branch->SetAddress(&fRecpoints);
b0e4d1e1 165
166 TBranch *branch1 = fTreeR->GetBranch("PMDRechit");
167 if (!branch1)
168 {
169 AliError("PMDRechit branch not found");
170 return;
171 }
172 branch1->SetAddress(&fRechits);
7dcaf913 173
26f1ae17 174 Int_t nmodules = (Int_t) branch->GetEntries();
175
ecee2a1a 176 AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
7dcaf913 177 for (Int_t imodule = 0; imodule < nmodules; imodule++)
178 {
26f1ae17 179 branch->GetEntry(imodule);
7dcaf913 180 Int_t nentries = fRecpoints->GetLast();
ecee2a1a 181 AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
182 ,nentries));
b0e4d1e1 183
184 Int_t ncrhit = 0;
185
7dcaf913 186 for(Int_t ient = 0; ient < nentries+1; ient++)
187 {
188 fPMDrecpoint = (AliPMDrecpoint1*)fRecpoints->UncheckedAt(ient);
189 idet = fPMDrecpoint->GetDetector();
190 ismn = fPMDrecpoint->GetSMNumber();
191 clusdata[0] = fPMDrecpoint->GetClusX();
192 clusdata[1] = fPMDrecpoint->GetClusY();
193 clusdata[2] = fPMDrecpoint->GetClusADC();
194 clusdata[3] = fPMDrecpoint->GetClusCells();
26f1ae17 195 clusdata[4] = fPMDrecpoint->GetClusSigmaX();
196 clusdata[5] = fPMDrecpoint->GetClusSigmaY();
920e13db 197
198 if (clusdata[4] != -99. && clusdata[5] != -99.)
199 {
b0e4d1e1 200 // extract the associated cell information
201 branch1->GetEntry(ncrhit);
202 Int_t nenbr1 = fRechits->GetLast() + 1;
203
204 irow = new Int_t[nenbr1];
205 icol = new Int_t[nenbr1];
206 itra = new Int_t[nenbr1];
207 ipid = new Int_t[nenbr1];
208 cadc = new Float_t[nenbr1];
209
210 for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
211 {
212 rechit = (AliPMDrechit*)fRechits->UncheckedAt(ient1);
213 //irow[ient1] = rechit->GetCellX();
214 //icol[ient1] = rechit->GetCellY();
215 itra[ient1] = rechit->GetCellTrack();
216 ipid[ient1] = rechit->GetCellPid();
217 cadc[ient1] = rechit->GetCellAdc();
218 }
219 AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
220 trackno, trackpid);
221
222 delete [] irow;
223 delete [] icol;
224 delete [] itra;
225 delete [] ipid;
226 delete [] cadc;
227
228 fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
920e13db 229 fPMDcontin->Add(fPMDclin);
b0e4d1e1 230
231 ncrhit++;
920e13db 232 }
7dcaf913 233 }
234 }
235
26f1ae17 236 AliPMDDiscriminator *pmddiscriminator = new AliPMDEmpDiscriminator();
237 pmddiscriminator->Discrimination(fPMDcontin,fPMDcontout);
7dcaf913 238
50555ba1 239 const Float_t kzpos = 361.5; // middle of the PMD
01c4d84a 240
b0e4d1e1 241 Int_t det,smn,trno,trpid,mstat;
7dcaf913 242 Float_t xpos,ypos;
b0e4d1e1 243 Float_t adc, ncell, radx, rady;
01c4d84a 244 Float_t xglobal = 0., yglobal = 0., zglobal = 0;
7dcaf913 245 Float_t pid;
246
7dcaf913 247
248 Int_t nentries2 = fPMDcontout->GetEntries();
ecee2a1a 249 AliDebug(1,Form("Number of clusters coming after discrimination = %d"
250 ,nentries2));
7dcaf913 251 for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
252 {
253 fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
254
255 det = fPMDclout->GetDetector();
256 smn = fPMDclout->GetSMN();
b0e4d1e1 257 trno = fPMDclout->GetClusTrackNo();
258 trpid = fPMDclout->GetClusTrackPid();
259 mstat = fPMDclout->GetClusMatching();
7dcaf913 260 xpos = fPMDclout->GetClusX();
261 ypos = fPMDclout->GetClusY();
262 adc = fPMDclout->GetClusADC();
263 ncell = fPMDclout->GetClusCells();
b0e4d1e1 264 radx = fPMDclout->GetClusSigmaX();
265 rady = fPMDclout->GetClusSigmaY();
7dcaf913 266 pid = fPMDclout->GetClusPID();
267
268 //
7dcaf913 269 /**********************************************************************
270 * det : Detector, 0: PRE & 1:CPV *
01c4d84a 271 * smn : Serial Module Number 0 to 23 for each plane *
7dcaf913 272 * xpos : x-position of the cluster *
273 * ypos : y-position of the cluster *
274 * THESE xpos & ypos are not the true xpos and ypos *
275 * for some of the unit modules. They are rotated. *
276 * adc : ADC contained in the cluster *
277 * ncell : Number of cells contained in the cluster *
278 * rad : radius of the cluster (1d fit) *
7dcaf913 279 **********************************************************************/
280 //
01c4d84a 281
282 fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal);
26f1ae17 283
284 if (det == 0)
285 {
50555ba1 286 zglobal = kzpos + 1.6; // PREshower plane
26f1ae17 287 }
288 else if (det == 1)
289 {
50555ba1 290 zglobal = kzpos - 1.7; // CPV plane
26f1ae17 291 }
292
7dcaf913 293 // Fill ESD
294
295 AliESDPmdTrack *esdpmdtr = new AliESDPmdTrack();
296
297 esdpmdtr->SetDetector(det);
c6f4d28b 298 esdpmdtr->SetSmn(smn);
299 esdpmdtr->SetClusterTrackNo(trno);
300 esdpmdtr->SetClusterTrackPid(trpid);
301 esdpmdtr->SetClusterMatching(mstat);
302
26f1ae17 303 esdpmdtr->SetClusterX(xglobal);
304 esdpmdtr->SetClusterY(yglobal);
305 esdpmdtr->SetClusterZ(zglobal);
7dcaf913 306 esdpmdtr->SetClusterADC(adc);
26f1ae17 307 esdpmdtr->SetClusterCells(ncell);
7dcaf913 308 esdpmdtr->SetClusterPID(pid);
c6f4d28b 309 esdpmdtr->SetClusterSigmaX(radx);
310 esdpmdtr->SetClusterSigmaY(rady);
7dcaf913 311
312 event->AddPmdTrack(esdpmdtr);
313 }
a1e9ac09 314
315 fPMDcontin->Delete();
316 fPMDcontout->Delete();
317
7dcaf913 318}
319//--------------------------------------------------------------------//
b0e4d1e1 320void AliPMDtracker::AssignTrPidToCluster(Int_t nentry, Int_t *itra,
321 Int_t *ipid, Float_t *cadc,
322 Int_t &trackno, Int_t &trackpid)
323{
324 // assign the track number and the corresponding pid to a cluster
325 // split cluster part will be done at the time of calculating eff/pur
326
327 Int_t *phentry = new Int_t [nentry];
328 Int_t *trenergy = 0x0;
329 Int_t *sortcoord = 0x0;
330
331 Int_t ngtrack = 0;
332 for (Int_t i = 0; i < nentry; i++)
333 {
334 phentry[i] = -1;
335 if (ipid[i] == 22)
336 {
337 phentry[ngtrack] = i;
338 ngtrack++;
339 }
340 }
341
342 if (ngtrack == 0)
343 {
344 // hadron track
345 // no need of track number, set to -1
346 trackpid = 8;
347 trackno = -1;
348 }
349 else if (ngtrack == 1)
350 {
351 // only one photon track
352 // track number set to photon track
353 trackpid = 1;
354 trackno = itra[phentry[0]];
355 }
356 else if (ngtrack > 1)
357 {
358 // more than one photon track
359
360 trenergy = new Int_t [ngtrack];
361 sortcoord = new Int_t [ngtrack];
362 for (Int_t i = 0; i < ngtrack; i++)
363 {
364 trenergy[i] = 0.;
365 for (Int_t j = 0; j < nentry; j++)
366 {
367 if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
368 {
369 trenergy[i] += cadc[j];
370 }
371 }
372 }
373
374 Bool_t jsort = true;
375 TMath::Sort(ngtrack,trenergy,sortcoord,jsort);
376
377 Int_t gtr = sortcoord[0];
378 trackno = itra[phentry[gtr]]; // highest adc track
379 trackpid = 1;
380
381 delete [] trenergy;
382 delete [] sortcoord;
383
384 } // end of ngtrack > 1
385
386
387
388}
389//--------------------------------------------------------------------//
7dcaf913 390void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
391{
392 fXvertex = vtx[0];
393 fYvertex = vtx[1];
394 fZvertex = vtx[2];
395 fSigmaX = evtx[0];
396 fSigmaY = evtx[1];
397 fSigmaZ = evtx[2];
398}
399//--------------------------------------------------------------------//
7dcaf913 400void AliPMDtracker::ResetClusters()
401{
402 if (fRecpoints) fRecpoints->Clear();
403}
404//--------------------------------------------------------------------//