Futher development of macros for AliEve preparations for PbPb
[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
be8b7039 145 Int_t idet = 0;
146 Int_t ismn = 0;
147 Int_t trackno = 1, trackpid = 0;
148 Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
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);
d270ca46 173
174 Int_t ncrhit = 0;
26f1ae17 175 Int_t nmodules = (Int_t) branch->GetEntries();
176
ecee2a1a 177 AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
7dcaf913 178 for (Int_t imodule = 0; imodule < nmodules; imodule++)
179 {
26f1ae17 180 branch->GetEntry(imodule);
7dcaf913 181 Int_t nentries = fRecpoints->GetLast();
ecee2a1a 182 AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
183 ,nentries));
b0e4d1e1 184
7dcaf913 185 for(Int_t ient = 0; ient < nentries+1; ient++)
186 {
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();
26f1ae17 194 clusdata[4] = fPMDrecpoint->GetClusSigmaX();
195 clusdata[5] = fPMDrecpoint->GetClusSigmaY();
4b1e197e 196
e6ba3040 197 if (clusdata[4] >= 0. && clusdata[5] >= 0.)
920e13db 198 {
b0e4d1e1 199 // extract the associated cell information
200 branch1->GetEntry(ncrhit);
201 Int_t nenbr1 = fRechits->GetLast() + 1;
202
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];
4b1e197e 208
b0e4d1e1 209 for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
210 {
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();
217 }
4b1e197e 218 if (idet == 0)
219 {
220 AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
221 trackno, trackpid);
222 }
223 else if (idet == 1)
224 {
225 trackno = itra[0];
226 trackpid = ipid[0];
227 }
b0e4d1e1 228
229 delete [] irow;
230 delete [] icol;
231 delete [] itra;
232 delete [] ipid;
233 delete [] cadc;
234
235 fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
920e13db 236 fPMDcontin->Add(fPMDclin);
b0e4d1e1 237
238 ncrhit++;
920e13db 239 }
7dcaf913 240 }
241 }
242
69c38808 243 AliPMDEmpDiscriminator pmddiscriminator;
244 pmddiscriminator.Discrimination(fPMDcontin,fPMDcontout);
7dcaf913 245
50555ba1 246 const Float_t kzpos = 361.5; // middle of the PMD
01c4d84a 247
e38d8aee 248 Int_t ix = -1, iy = -1;
be8b7039 249 Int_t det = 0, smn = 0, trno = 1, trpid = 0, mstat = 0;
250 Float_t xpos = 0., ypos = 0.;
251 Float_t adc = 0., ncell = 0., radx = 0., rady = 0.;
01c4d84a 252 Float_t xglobal = 0., yglobal = 0., zglobal = 0;
be8b7039 253 Float_t pid = 0.;
7dcaf913 254
dc461f61 255 fPMDutil->ApplyAlignment();
7dcaf913 256
257 Int_t nentries2 = fPMDcontout->GetEntries();
ecee2a1a 258 AliDebug(1,Form("Number of clusters coming after discrimination = %d"
259 ,nentries2));
7dcaf913 260 for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
261 {
262 fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
263
264 det = fPMDclout->GetDetector();
265 smn = fPMDclout->GetSMN();
b0e4d1e1 266 trno = fPMDclout->GetClusTrackNo();
267 trpid = fPMDclout->GetClusTrackPid();
268 mstat = fPMDclout->GetClusMatching();
7dcaf913 269 xpos = fPMDclout->GetClusX();
270 ypos = fPMDclout->GetClusY();
271 adc = fPMDclout->GetClusADC();
272 ncell = fPMDclout->GetClusCells();
b0e4d1e1 273 radx = fPMDclout->GetClusSigmaX();
6db43724 274 // Here in the variable "rady" we are keeping the row and col
275 // of the single isolated cluster having ncell = 1 for offline
276 // calibration
277
278 if ((radx > 999. && radx < 1000.) && ncell == 1)
279 {
280 if (smn < 12)
281 {
282 ix = (Int_t) (ypos +0.5);
283 iy = (Int_t) xpos;
284 }
285 else if (smn > 12 && smn < 24)
286 {
287 ix = (Int_t) xpos;
288 iy = (Int_t) (ypos +0.5);
289 }
290 rady = (Float_t) (ix*100 + iy);
291 }
292 else
293 {
294 rady = fPMDclout->GetClusSigmaY();
295 }
7dcaf913 296 pid = fPMDclout->GetClusPID();
297
298 //
7dcaf913 299 /**********************************************************************
300 * det : Detector, 0: PRE & 1:CPV *
01c4d84a 301 * smn : Serial Module Number 0 to 23 for each plane *
7dcaf913 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) *
7dcaf913 309 **********************************************************************/
310 //
01c4d84a 311
26f1ae17 312
313 if (det == 0)
314 {
dc461f61 315 zglobal = kzpos + 1.65; // PREshower plane
26f1ae17 316 }
317 else if (det == 1)
318 {
dc461f61 319 zglobal = kzpos - 1.65; // CPV plane
26f1ae17 320 }
321
dc461f61 322 fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal,zglobal);
323
7dcaf913 324 // Fill ESD
325
326 AliESDPmdTrack *esdpmdtr = new AliESDPmdTrack();
327
328 esdpmdtr->SetDetector(det);
c6f4d28b 329 esdpmdtr->SetSmn(smn);
330 esdpmdtr->SetClusterTrackNo(trno);
331 esdpmdtr->SetClusterTrackPid(trpid);
332 esdpmdtr->SetClusterMatching(mstat);
333
26f1ae17 334 esdpmdtr->SetClusterX(xglobal);
335 esdpmdtr->SetClusterY(yglobal);
336 esdpmdtr->SetClusterZ(zglobal);
7dcaf913 337 esdpmdtr->SetClusterADC(adc);
26f1ae17 338 esdpmdtr->SetClusterCells(ncell);
7dcaf913 339 esdpmdtr->SetClusterPID(pid);
c6f4d28b 340 esdpmdtr->SetClusterSigmaX(radx);
341 esdpmdtr->SetClusterSigmaY(rady);
7dcaf913 342
343 event->AddPmdTrack(esdpmdtr);
69c38808 344 delete esdpmdtr;
7dcaf913 345 }
a1e9ac09 346
347 fPMDcontin->Delete();
348 fPMDcontout->Delete();
349
7dcaf913 350}
351//--------------------------------------------------------------------//
b0e4d1e1 352void AliPMDtracker::AssignTrPidToCluster(Int_t nentry, Int_t *itra,
353 Int_t *ipid, Float_t *cadc,
354 Int_t &trackno, Int_t &trackpid)
355{
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
358
359 Int_t *phentry = new Int_t [nentry];
a6621fb4 360 Int_t *hadentry = new Int_t [nentry];
b0e4d1e1 361 Int_t *trenergy = 0x0;
a6621fb4 362 Int_t *trpid = 0x0;
b0e4d1e1 363 Int_t *sortcoord = 0x0;
364
365 Int_t ngtrack = 0;
a6621fb4 366 Int_t nhtrack = 0;
b0e4d1e1 367 for (Int_t i = 0; i < nentry; i++)
368 {
369 phentry[i] = -1;
a6621fb4 370 hadentry[i] = -1;
371
b0e4d1e1 372 if (ipid[i] == 22)
373 {
374 phentry[ngtrack] = i;
375 ngtrack++;
376 }
a6621fb4 377 else if (ipid[i] != 22)
378 {
379 hadentry[nhtrack] = i;
380 nhtrack++;
381 }
b0e4d1e1 382 }
4b1e197e 383
a6621fb4 384 Int_t nghadtrack = ngtrack + nhtrack;
385
b0e4d1e1 386 if (ngtrack == 0)
387 {
388 // hadron track
389 // no need of track number, set to -1
390 trackpid = 8;
391 trackno = -1;
392 }
a6621fb4 393 else if (ngtrack >= 1)
b0e4d1e1 394 {
a6621fb4 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
398
399 trenergy = new Int_t [nghadtrack];
400 trpid = new Int_t [nghadtrack];
401 sortcoord = new Int_t [nghadtrack];
b0e4d1e1 402 for (Int_t i = 0; i < ngtrack; i++)
403 {
404 trenergy[i] = 0.;
a6621fb4 405 trpid[i] = -1;
b0e4d1e1 406 for (Int_t j = 0; j < nentry; j++)
407 {
408 if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
409 {
410 trenergy[i] += cadc[j];
a6621fb4 411 trpid[i] = 22;
412 }
413 }
414 }
415 for (Int_t i = ngtrack; i < nghadtrack; i++)
416 {
417 trenergy[i] = 0.;
418 trpid[i] = -1;
419 for (Int_t j = 0; j < nentry; j++)
420 {
421 if (ipid[j] != 22 && itra[j] == itra[hadentry[i-ngtrack]])
422 {
423 trenergy[i] += cadc[j];
424 trpid[i] = ipid[j];
b0e4d1e1 425 }
426 }
427 }
4b1e197e 428
b0e4d1e1 429 Bool_t jsort = true;
a6621fb4 430 TMath::Sort(nghadtrack,trenergy,sortcoord,jsort);
4b1e197e 431
b0e4d1e1 432 Int_t gtr = sortcoord[0];
a6621fb4 433 if (trpid[gtr] == 22)
434 {
435 trackpid = 22;
436 trackno = itra[phentry[gtr]]; // highest adc track
437 }
438 else
439 {
440 trackpid = 8;
441 trackno = -1;
442 }
4b1e197e 443
b0e4d1e1 444 delete [] trenergy;
a6621fb4 445 delete [] trpid;
b0e4d1e1 446 delete [] sortcoord;
4b1e197e 447
a6621fb4 448 } // end of ngtrack >= 1
4b1e197e 449
b0e4d1e1 450}
451//--------------------------------------------------------------------//
7dcaf913 452void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
453{
454 fXvertex = vtx[0];
455 fYvertex = vtx[1];
456 fZvertex = vtx[2];
457 fSigmaX = evtx[0];
458 fSigmaY = evtx[1];
459 fSigmaZ = evtx[2];
460}
461//--------------------------------------------------------------------//
7dcaf913 462void AliPMDtracker::ResetClusters()
463{
464 if (fRecpoints) fRecpoints->Clear();
465}
466//--------------------------------------------------------------------//