Update in gain cdb macro: sjena
[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
32b0cad8 35#include <TGeoMatrix.h>
36
37#include "AliGeomManager.h"
38
7dcaf913 39#include "AliPMDcluster.h"
40#include "AliPMDclupid.h"
41#include "AliPMDrecpoint1.h"
b0e4d1e1 42#include "AliPMDrecdata.h"
43#include "AliPMDrechit.h"
7dcaf913 44#include "AliPMDUtility.h"
45#include "AliPMDDiscriminator.h"
26f1ae17 46#include "AliPMDEmpDiscriminator.h"
7dcaf913 47#include "AliPMDtracker.h"
48
49#include "AliESDPmdTrack.h"
af885e0f 50#include "AliESDEvent.h"
ecee2a1a 51#include "AliLog.h"
7dcaf913 52
53ClassImp(AliPMDtracker)
54
55AliPMDtracker::AliPMDtracker():
56 fTreeR(0),
722ccc67 57 fRecpoints(new TClonesArray("AliPMDrecpoint1", 10)),
b0e4d1e1 58 fRechits(new TClonesArray("AliPMDrechit", 10)),
7dcaf913 59 fPMDcontin(new TObjArray()),
60 fPMDcontout(new TObjArray()),
7dcaf913 61 fPMDutil(new AliPMDUtility()),
62 fPMDrecpoint(0),
63 fPMDclin(0),
64 fPMDclout(0),
7dcaf913 65 fXvertex(0.),
66 fYvertex(0.),
67 fZvertex(0.),
68 fSigmaX(0.),
69 fSigmaY(0.),
70 fSigmaZ(0.)
71{
72 //
73 // Default Constructor
74 //
75}
76//--------------------------------------------------------------------//
a48edddd 77AliPMDtracker:: AliPMDtracker(const AliPMDtracker & /* tracker */):
2332574a 78 TObject(/* tracker */),
79 fTreeR(0),
80 fRecpoints(NULL),
b0e4d1e1 81 fRechits(NULL),
2332574a 82 fPMDcontin(NULL),
83 fPMDcontout(NULL),
84 fPMDutil(NULL),
85 fPMDrecpoint(0),
86 fPMDclin(0),
87 fPMDclout(0),
88 fXvertex(0.),
89 fYvertex(0.),
90 fZvertex(0.),
91 fSigmaX(0.),
92 fSigmaY(0.),
93 fSigmaZ(0.)
a48edddd 94{
95 // copy constructor
96 AliError("Copy constructor not allowed");
97}
98
99//--------------------------------------------------------------------//
100AliPMDtracker& AliPMDtracker::operator=(const AliPMDtracker & /* tracker */)
101{
102 // assignment operator
103 AliError("Assignment operator not allowed");
104 return *this;
105}
106
107//--------------------------------------------------------------------//
7dcaf913 108AliPMDtracker::~AliPMDtracker()
109{
110 // Destructor
111 if (fRecpoints)
112 {
722ccc67 113 fRecpoints->Clear();
7dcaf913 114 }
b0e4d1e1 115 if (fRechits)
116 {
117 fRechits->Clear();
118 }
119
7dcaf913 120 if (fPMDcontin)
121 {
122 fPMDcontin->Delete();
123 delete fPMDcontin;
124 fPMDcontin=0;
a1e9ac09 125
7dcaf913 126 }
127 if (fPMDcontout)
722ccc67 128 {
7dcaf913 129 fPMDcontout->Delete();
130 delete fPMDcontout;
131 fPMDcontout=0;
a1e9ac09 132
7dcaf913 133 }
316c6cd9 134 delete fPMDutil;
7dcaf913 135}
136//--------------------------------------------------------------------//
137void AliPMDtracker::LoadClusters(TTree *treein)
138{
139 // Load the Reconstructed tree
140 fTreeR = treein;
141}
142//--------------------------------------------------------------------//
af885e0f 143void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
7dcaf913 144{
145 // Converts digits to recpoints after running clustering
146 // algorithm on CPV plane and PREshower plane
147 //
148
be8b7039 149 Int_t idet = 0;
150 Int_t ismn = 0;
151 Int_t trackno = 1, trackpid = 0;
152 Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
b0e4d1e1 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;
7dcaf913 161
162 TBranch *branch = fTreeR->GetBranch("PMDRecpoint");
ecee2a1a 163 if (!branch)
164 {
165 AliError("PMDRecpoint branch not found");
166 return;
167 }
7dcaf913 168 branch->SetAddress(&fRecpoints);
b0e4d1e1 169
170 TBranch *branch1 = fTreeR->GetBranch("PMDRechit");
171 if (!branch1)
172 {
173 AliError("PMDRechit branch not found");
174 return;
175 }
176 branch1->SetAddress(&fRechits);
d270ca46 177
178 Int_t ncrhit = 0;
26f1ae17 179 Int_t nmodules = (Int_t) branch->GetEntries();
180
ecee2a1a 181 AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
7dcaf913 182 for (Int_t imodule = 0; imodule < nmodules; imodule++)
183 {
26f1ae17 184 branch->GetEntry(imodule);
7dcaf913 185 Int_t nentries = fRecpoints->GetLast();
ecee2a1a 186 AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
187 ,nentries));
b0e4d1e1 188
7dcaf913 189 for(Int_t ient = 0; ient < nentries+1; ient++)
190 {
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();
26f1ae17 198 clusdata[4] = fPMDrecpoint->GetClusSigmaX();
199 clusdata[5] = fPMDrecpoint->GetClusSigmaY();
4b1e197e 200
e6ba3040 201 if (clusdata[4] >= 0. && clusdata[5] >= 0.)
920e13db 202 {
b0e4d1e1 203 // extract the associated cell information
204 branch1->GetEntry(ncrhit);
205 Int_t nenbr1 = fRechits->GetLast() + 1;
206
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];
4b1e197e 212
b0e4d1e1 213 for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
214 {
db06ef51 215 irow[ient1] = -99;
216 icol[ient1] = -99;
217 itra[ient1] = -99;
218 ipid[ient1] = -99;
219 cadc[ient1] = 0.;
220 }
221 for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
222 {
b0e4d1e1 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();
229 }
4b1e197e 230 if (idet == 0)
231 {
232 AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
233 trackno, trackpid);
234 }
235 else if (idet == 1)
236 {
237 trackno = itra[0];
238 trackpid = ipid[0];
239 }
b0e4d1e1 240
241 delete [] irow;
242 delete [] icol;
243 delete [] itra;
244 delete [] ipid;
245 delete [] cadc;
246
247 fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
920e13db 248 fPMDcontin->Add(fPMDclin);
b0e4d1e1 249
250 ncrhit++;
920e13db 251 }
7dcaf913 252 }
253 }
254
69c38808 255 AliPMDEmpDiscriminator pmddiscriminator;
256 pmddiscriminator.Discrimination(fPMDcontin,fPMDcontout);
7dcaf913 257
32b0cad8 258 // alignment implemention
259
260 Double_t sectr[4][3] = { {0.,0.,0.},{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
261 TString snsector="PMD/Sector";
262 TString symname;
263 TGeoHMatrix gpmdor;
264
265 for(Int_t isector=1; isector<=4; isector++)
266 {
267 symname = snsector;
268 symname += isector;
269 TGeoHMatrix *gpmdal = AliGeomManager::GetMatrix(symname);
270 Double_t *tral = gpmdal->GetTranslation();
271
272 AliGeomManager::GetOrigGlobalMatrix(symname, gpmdor);
273 Double_t *tror = gpmdor.GetTranslation();
274
275 for(Int_t ixyz=0; ixyz<3; ixyz++)
276 {
277 sectr[isector-1][ixyz] = tral[ixyz] - tror[ixyz];
278 }
279 }
280
50555ba1 281 const Float_t kzpos = 361.5; // middle of the PMD
01c4d84a 282
e38d8aee 283 Int_t ix = -1, iy = -1;
be8b7039 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.;
01c4d84a 287 Float_t xglobal = 0., yglobal = 0., zglobal = 0;
be8b7039 288 Float_t pid = 0.;
7dcaf913 289
32b0cad8 290 fPMDutil->ApplyAlignment(sectr);
7dcaf913 291
292 Int_t nentries2 = fPMDcontout->GetEntries();
ecee2a1a 293 AliDebug(1,Form("Number of clusters coming after discrimination = %d"
294 ,nentries2));
7dcaf913 295 for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
296 {
297 fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
298
299 det = fPMDclout->GetDetector();
300 smn = fPMDclout->GetSMN();
b0e4d1e1 301 trno = fPMDclout->GetClusTrackNo();
302 trpid = fPMDclout->GetClusTrackPid();
303 mstat = fPMDclout->GetClusMatching();
7dcaf913 304 xpos = fPMDclout->GetClusX();
305 ypos = fPMDclout->GetClusY();
306 adc = fPMDclout->GetClusADC();
307 ncell = fPMDclout->GetClusCells();
b0e4d1e1 308 radx = fPMDclout->GetClusSigmaX();
6db43724 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 if ((radx > 999. && radx < 1000.) && ncell == 1)
314 {
315 if (smn < 12)
316 {
317 ix = (Int_t) (ypos +0.5);
318 iy = (Int_t) xpos;
319 }
320 else if (smn > 12 && smn < 24)
321 {
322 ix = (Int_t) xpos;
323 iy = (Int_t) (ypos +0.5);
324 }
325 rady = (Float_t) (ix*100 + iy);
326 }
327 else
328 {
329 rady = fPMDclout->GetClusSigmaY();
330 }
7dcaf913 331 pid = fPMDclout->GetClusPID();
332
333 //
7dcaf913 334 /**********************************************************************
335 * det : Detector, 0: PRE & 1:CPV *
01c4d84a 336 * smn : Serial Module Number 0 to 23 for each plane *
7dcaf913 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) *
7dcaf913 344 **********************************************************************/
345 //
01c4d84a 346
26f1ae17 347
348 if (det == 0)
349 {
dc461f61 350 zglobal = kzpos + 1.65; // PREshower plane
26f1ae17 351 }
352 else if (det == 1)
353 {
dc461f61 354 zglobal = kzpos - 1.65; // CPV plane
26f1ae17 355 }
356
dc461f61 357 fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal,zglobal);
358
7dcaf913 359 // Fill ESD
360
361 AliESDPmdTrack *esdpmdtr = new AliESDPmdTrack();
362
363 esdpmdtr->SetDetector(det);
c6f4d28b 364 esdpmdtr->SetSmn(smn);
365 esdpmdtr->SetClusterTrackNo(trno);
366 esdpmdtr->SetClusterTrackPid(trpid);
367 esdpmdtr->SetClusterMatching(mstat);
368
26f1ae17 369 esdpmdtr->SetClusterX(xglobal);
370 esdpmdtr->SetClusterY(yglobal);
371 esdpmdtr->SetClusterZ(zglobal);
7dcaf913 372 esdpmdtr->SetClusterADC(adc);
26f1ae17 373 esdpmdtr->SetClusterCells(ncell);
7dcaf913 374 esdpmdtr->SetClusterPID(pid);
c6f4d28b 375 esdpmdtr->SetClusterSigmaX(radx);
376 esdpmdtr->SetClusterSigmaY(rady);
7dcaf913 377
378 event->AddPmdTrack(esdpmdtr);
69c38808 379 delete esdpmdtr;
7dcaf913 380 }
a1e9ac09 381
382 fPMDcontin->Delete();
383 fPMDcontout->Delete();
384
7dcaf913 385}
386//--------------------------------------------------------------------//
b0e4d1e1 387void 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 Int_t *phentry = new Int_t [nentry];
a6621fb4 395 Int_t *hadentry = new Int_t [nentry];
a6621fb4 396 Int_t *trpid = 0x0;
b0e4d1e1 397 Int_t *sortcoord = 0x0;
b00d6574 398 Float_t *trenergy = 0x0;
b0e4d1e1 399
400 Int_t ngtrack = 0;
a6621fb4 401 Int_t nhtrack = 0;
b0e4d1e1 402 for (Int_t i = 0; i < nentry; i++)
403 {
404 phentry[i] = -1;
a6621fb4 405 hadentry[i] = -1;
406
b0e4d1e1 407 if (ipid[i] == 22)
408 {
409 phentry[ngtrack] = i;
410 ngtrack++;
411 }
a6621fb4 412 else if (ipid[i] != 22)
413 {
414 hadentry[nhtrack] = i;
415 nhtrack++;
416 }
b0e4d1e1 417 }
4b1e197e 418
a6621fb4 419 Int_t nghadtrack = ngtrack + nhtrack;
420
b0e4d1e1 421 if (ngtrack == 0)
422 {
423 // hadron track
424 // no need of track number, set to -1
425 trackpid = 8;
426 trackno = -1;
427 }
a6621fb4 428 else if (ngtrack >= 1)
b0e4d1e1 429 {
a6621fb4 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
b00d6574 434 trenergy = new Float_t [nghadtrack];
a6621fb4 435 trpid = new Int_t [nghadtrack];
47a8a50d 436 sortcoord = new Int_t [2*nghadtrack];
b0e4d1e1 437 for (Int_t i = 0; i < ngtrack; i++)
438 {
439 trenergy[i] = 0.;
a6621fb4 440 trpid[i] = -1;
b0e4d1e1 441 for (Int_t j = 0; j < nentry; j++)
442 {
443 if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
444 {
445 trenergy[i] += cadc[j];
a6621fb4 446 trpid[i] = 22;
447 }
448 }
449 }
450 for (Int_t i = ngtrack; i < nghadtrack; i++)
451 {
452 trenergy[i] = 0.;
453 trpid[i] = -1;
454 for (Int_t j = 0; j < nentry; j++)
455 {
456 if (ipid[j] != 22 && itra[j] == itra[hadentry[i-ngtrack]])
457 {
458 trenergy[i] += cadc[j];
459 trpid[i] = ipid[j];
b0e4d1e1 460 }
461 }
462 }
4b1e197e 463
b0e4d1e1 464 Bool_t jsort = true;
a6621fb4 465 TMath::Sort(nghadtrack,trenergy,sortcoord,jsort);
4b1e197e 466
b0e4d1e1 467 Int_t gtr = sortcoord[0];
a6621fb4 468 if (trpid[gtr] == 22)
469 {
470 trackpid = 22;
471 trackno = itra[phentry[gtr]]; // highest adc track
472 }
473 else
474 {
475 trackpid = 8;
476 trackno = -1;
477 }
4b1e197e 478
b0e4d1e1 479 delete [] trenergy;
a6621fb4 480 delete [] trpid;
b0e4d1e1 481 delete [] sortcoord;
4b1e197e 482
a6621fb4 483 } // end of ngtrack >= 1
db06ef51 484
485 delete [] phentry;
486 delete [] hadentry;
4b1e197e 487
b0e4d1e1 488}
489//--------------------------------------------------------------------//
7dcaf913 490void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
491{
492 fXvertex = vtx[0];
493 fYvertex = vtx[1];
494 fZvertex = vtx[2];
495 fSigmaX = evtx[0];
496 fSigmaY = evtx[1];
497 fSigmaZ = evtx[2];
498}
499//--------------------------------------------------------------------//
7dcaf913 500void AliPMDtracker::ResetClusters()
501{
502 if (fRecpoints) fRecpoints->Clear();
503}
504//--------------------------------------------------------------------//