Possibility to change sign of Pz event by event, needed for pA, Ap simulation in
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizerV1.cxx
CommitLineData
f7336fa3 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
88cb7938 16/* $Id$ */
f7336fa3 17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// TRD cluster finder for the slow simulator.
21// //
22///////////////////////////////////////////////////////////////////////////////
23
24#include <TF1.h>
94de3818 25#include <TTree.h>
793ff80c 26#include <TH1.h>
a819a5f7 27#include <TFile.h>
f7336fa3 28
793ff80c 29#include "AliRun.h"
88cb7938 30#include "AliRunLoader.h"
31#include "AliLoader.h"
793ff80c 32
f7336fa3 33#include "AliTRDclusterizerV1.h"
34#include "AliTRDmatrix.h"
35#include "AliTRDgeometry.h"
6f1e466d 36#include "AliTRDdataArrayF.h"
793ff80c 37#include "AliTRDdataArrayI.h"
38#include "AliTRDdigitsManager.h"
17b26de4 39#include "AliTRDparameter.h"
f7336fa3 40
41ClassImp(AliTRDclusterizerV1)
42
43//_____________________________________________________________________________
44AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
45{
46 //
47 // AliTRDclusterizerV1 default constructor
48 //
49
17b26de4 50 fDigitsManager = 0;
db30bf0f 51
f7336fa3 52}
53
54//_____________________________________________________________________________
55AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
56 :AliTRDclusterizer(name,title)
57{
58 //
59 // AliTRDclusterizerV1 default constructor
60 //
61
6f1e466d 62 fDigitsManager = new AliTRDdigitsManager();
17b26de4 63 fDigitsManager->CreateArrays();
f7336fa3 64
65}
66
8230f242 67//_____________________________________________________________________________
dd9a6ee3 68AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
73ae7b59 69:AliTRDclusterizer(c)
8230f242 70{
71 //
72 // AliTRDclusterizerV1 copy constructor
73 //
74
dd9a6ee3 75 ((AliTRDclusterizerV1 &) c).Copy(*this);
8230f242 76
77}
78
f7336fa3 79//_____________________________________________________________________________
80AliTRDclusterizerV1::~AliTRDclusterizerV1()
81{
8230f242 82 //
83 // AliTRDclusterizerV1 destructor
84 //
f7336fa3 85
6f1e466d 86 if (fDigitsManager) {
87 delete fDigitsManager;
abaf1f1d 88 fDigitsManager = NULL;
f7336fa3 89 }
90
91}
92
dd9a6ee3 93//_____________________________________________________________________________
94AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
95{
96 //
97 // Assignment operator
98 //
99
100 if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
101 return *this;
102
103}
104
8230f242 105//_____________________________________________________________________________
e0d47c25 106void AliTRDclusterizerV1::Copy(TObject &c) const
8230f242 107{
108 //
109 // Copy function
110 //
111
17b26de4 112 ((AliTRDclusterizerV1 &) c).fDigitsManager = 0;
8230f242 113
114 AliTRDclusterizer::Copy(c);
115
116}
117
f7336fa3 118//_____________________________________________________________________________
119Bool_t AliTRDclusterizerV1::ReadDigits()
120{
121 //
122 // Reads the digits arrays from the input aliroot file
123 //
124
88cb7938 125 if (!fRunLoader) {
17b26de4 126 printf("<AliTRDclusterizerV1::ReadDigits> ");
f7336fa3 127 printf("No input file open\n");
128 return kFALSE;
129 }
88cb7938 130 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
131 if (!loader->TreeD()) loader->LoadDigits();
abaf1f1d 132
f7336fa3 133 // Read in the digit arrays
88cb7938 134 return (fDigitsManager->ReadDigits(loader->TreeD()));
f7336fa3 135
136}
137
138//_____________________________________________________________________________
793ff80c 139Bool_t AliTRDclusterizerV1::MakeClusters()
f7336fa3 140{
141 //
142 // Generates the cluster.
143 //
144
145 Int_t row, col, time;
146
bdbb05bb 147 /*
3e1a3ad8 148 if (fTRD->IsVersion() != 1) {
17b26de4 149 printf("<AliTRDclusterizerV1::MakeCluster> ");
f7336fa3 150 printf("TRD must be version 1 (slow simulator).\n");
151 return kFALSE;
152 }
bdbb05bb 153 */
f7336fa3 154
155 // Get the geometry
bdbb05bb 156 AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader);
f7336fa3 157
17b26de4 158 // Create a default parameter class if none is defined
159 if (!fPar) {
160 fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
5443e65e 161 printf("<AliTRDclusterizerV1::MakeCluster> ");
162 printf("Create the default parameter object.\n");
17b26de4 163 }
164
ccb4315c 165 Float_t timeBinSize = fPar->GetDriftVelocity()
166 / fPar->GetSamplingFrequency();
a819a5f7 167 // Half of ampl.region
168 const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.;
169
17b26de4 170 Float_t omegaTau = fPar->GetOmegaTau();
47517f42 171 if (fVerbose > 0) {
17b26de4 172 printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42 173 printf("OmegaTau = %f \n",omegaTau);
17b26de4 174 printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42 175 printf("Start creating clusters.\n");
176 }
f7336fa3 177
8230f242 178 AliTRDdataArrayI *digits;
793ff80c 179 AliTRDdataArrayI *track0;
180 AliTRDdataArrayI *track1;
181 AliTRDdataArrayI *track2;
f7336fa3 182
3e1a3ad8 183 // Threshold value for the maximum
17b26de4 184 Int_t maxThresh = fPar->GetClusMaxThresh();
3e1a3ad8 185 // Threshold value for the digit signal
17b26de4 186 Int_t sigThresh = fPar->GetClusSigThresh();
f7336fa3 187
188 // Iteration limit for unfolding procedure
8230f242 189 const Float_t kEpsilon = 0.01;
f7336fa3 190
8230f242 191 const Int_t kNclus = 3;
192 const Int_t kNsig = 5;
3e1a3ad8 193 const Int_t kNtrack = 3 * kNclus;
194
db30bf0f 195 Int_t iType = 0;
196 Int_t iUnfold = 0;
197
198 Float_t ratioLeft = 1.0;
199 Float_t ratioRight = 1.0;
200
3e1a3ad8 201 Float_t padSignal[kNsig];
202 Float_t clusterSignal[kNclus];
203 Float_t clusterPads[kNclus];
204 Int_t clusterDigit[kNclus];
205 Int_t clusterTracks[kNtrack];
f7336fa3 206
207 Int_t chamBeg = 0;
793ff80c 208 Int_t chamEnd = AliTRDgeometry::Ncham();
f7336fa3 209 Int_t planBeg = 0;
793ff80c 210 Int_t planEnd = AliTRDgeometry::Nplan();
f7336fa3 211 Int_t sectBeg = 0;
793ff80c 212 Int_t sectEnd = AliTRDgeometry::Nsect();
f7336fa3 213
3e1a3ad8 214 // Start clustering in every chamber
f7336fa3 215 for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
216 for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
217 for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
218
8230f242 219 Int_t idet = geo->GetDetector(iplan,icham,isect);
f7336fa3 220
db30bf0f 221 Int_t nClusters = 0;
222 Int_t nClusters2pad = 0;
223 Int_t nClusters3pad = 0;
224 Int_t nClusters4pad = 0;
225 Int_t nClusters5pad = 0;
226 Int_t nClustersLarge = 0;
3e1a3ad8 227
47517f42 228 if (fVerbose > 0) {
17b26de4 229 printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42 230 printf("Analyzing chamber %d, plane %d, sector %d.\n"
231 ,icham,iplan,isect);
232 }
f7336fa3 233
5443e65e 234 Int_t nRowMax = fPar->GetRowMax(iplan,icham,isect);
235 Int_t nColMax = fPar->GetColMax(iplan);
236 Int_t nTimeBefore = fPar->GetTimeBefore();
237 Int_t nTimeTotal = fPar->GetTimeTotal();
238
239 Float_t row0 = fPar->GetRow0(iplan,icham,isect);
240 Float_t col0 = fPar->GetCol0(iplan);
241 Float_t rowSize = fPar->GetRowPadSize(iplan,icham,isect);
242 Float_t colSize = fPar->GetColPadSize(iplan);
f7336fa3 243
3e1a3ad8 244 // Get the digits
8230f242 245 digits = fDigitsManager->GetDigits(idet);
3e1a3ad8 246 digits->Expand();
793ff80c 247 track0 = fDigitsManager->GetDictionary(idet,0);
3e1a3ad8 248 track0->Expand();
793ff80c 249 track1 = fDigitsManager->GetDictionary(idet,1);
3e1a3ad8 250 track1->Expand();
793ff80c 251 track2 = fDigitsManager->GetDictionary(idet,2);
3e1a3ad8 252 track2->Expand();
253
254 // Loop through the chamber and find the maxima
255 for ( row = 0; row < nRowMax; row++) {
95d99e62 256 // for ( col = 2; col < nColMax; col++) {
257 for ( col = 4; col < nColMax-2; col++) {
3e1a3ad8 258 for (time = 0; time < nTimeTotal; time++) {
259
a819a5f7 260 Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col ,time));
261 Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
262 Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
3e1a3ad8 263
264 // Look for the maximum
db30bf0f 265 if (signalM >= maxThresh) {
266 if (((signalL >= sigThresh) &&
267 (signalL < signalM)) ||
268 ((signalR >= sigThresh) &&
269 (signalR < signalM))) {
3e1a3ad8 270 // Maximum found, mark the position by a negative signal
271 digits->SetDataUnchecked(row,col-1,time,-signalM);
272 }
273 }
274
275 }
276 }
277 }
278
279 // Now check the maxima and calculate the cluster position
280 for ( row = 0; row < nRowMax ; row++) {
db30bf0f 281 for (time = 0; time < nTimeTotal; time++) {
282 for ( col = 1; col < nColMax-1; col++) {
3e1a3ad8 283
284 // Maximum found ?
285 if (digits->GetDataUnchecked(row,col,time) < 0) {
f7336fa3 286
9d0b222b 287 Int_t iPad;
8230f242 288 for (iPad = 0; iPad < kNclus; iPad++) {
3e1a3ad8 289 Int_t iPadCol = col - 1 + iPad;
290 clusterSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
291 ,iPadCol
292 ,time));
293 clusterDigit[iPad] = digits->GetIndexUnchecked(row,iPadCol,time);
294 clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
295 clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
296 clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
f7336fa3 297 }
298
db30bf0f 299 // Count the number of pads in the cluster
300 Int_t nPadCount = 0;
301 Int_t ii = 0;
302 while (TMath::Abs(digits->GetDataUnchecked(row,col-ii ,time))
303 >= sigThresh) {
304 nPadCount++;
305 ii++;
306 if (col-ii < 0) break;
307 }
308 ii = 0;
309 while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time))
310 >= sigThresh) {
311 nPadCount++;
312 ii++;
313 if (col+ii+1 >= nColMax) break;
314 }
315
316 nClusters++;
317 switch (nPadCount) {
318 case 2:
319 iType = 0;
320 nClusters2pad++;
321 break;
322 case 3:
323 iType = 1;
324 nClusters3pad++;
325 break;
326 case 4:
327 iType = 2;
328 nClusters4pad++;
329 break;
330 case 5:
331 iType = 3;
332 nClusters5pad++;
333 break;
334 default:
335 iType = 4;
336 nClustersLarge++;
337 break;
338 };
339
340 // Don't analyze large clusters
341 //if (iType == 4) continue;
342
343 // Look for 5 pad cluster with minimum in the middle
344 Bool_t fivePadCluster = kFALSE;
3e1a3ad8 345 if (col < nColMax-3) {
346 if (digits->GetDataUnchecked(row,col+2,time) < 0) {
db30bf0f 347 fivePadCluster = kTRUE;
348 }
349 if ((fivePadCluster) && (col < nColMax-5)) {
350 if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) {
351 fivePadCluster = kFALSE;
352 }
353 }
354 if ((fivePadCluster) && (col > 1)) {
355 if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) {
356 fivePadCluster = kFALSE;
357 }
358 }
359 }
360
361 // 5 pad cluster
362 // Modify the signal of the overlapping pad for the left part
363 // of the cluster which remains from a previous unfolding
364 if (iUnfold) {
365 clusterSignal[0] *= ratioLeft;
366 iType = 3;
367 iUnfold = 0;
368 }
369
370 // Unfold the 5 pad cluster
371 if (fivePadCluster) {
372 for (iPad = 0; iPad < kNsig; iPad++) {
373 padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
374 ,col-1+iPad
375 ,time));
f7336fa3 376 }
db30bf0f 377 // Unfold the two maxima and set the signal on
378 // the overlapping pad to the ratio
17b26de4 379 ratioRight = Unfold(kEpsilon,iplan,padSignal);
db30bf0f 380 ratioLeft = 1.0 - ratioRight;
381 clusterSignal[2] *= ratioRight;
382 iType = 3;
383 iUnfold = 1;
f7336fa3 384 }
f7336fa3 385
3e1a3ad8 386 Float_t clusterCharge = clusterSignal[0]
387 + clusterSignal[1]
388 + clusterSignal[2];
389
db30bf0f 390 // The position of the cluster
3e1a3ad8 391 clusterPads[0] = row + 0.5;
3e1a3ad8 392 // Take the shift of the additional time bins into account
393 clusterPads[2] = time - nTimeBefore + 0.5;
394
17b26de4 395 if (fPar->LUTOn()) {
db30bf0f 396
397 // Calculate the position of the cluster by using the
398 // lookup table method
5443e65e 399 clusterPads[1] = col + 0.5
400 + fPar->LUTposition(iplan,clusterSignal[0]
17b26de4 401 ,clusterSignal[1]
402 ,clusterSignal[2]);
db30bf0f 403
404 }
405 else {
406
407 // Calculate the position of the cluster by using the
408 // center of gravity method
409 clusterPads[1] = col + 0.5
410 + (clusterSignal[2] - clusterSignal[0])
411 / clusterCharge;
412
413 }
414
5443e65e 415 Float_t q0 = clusterSignal[0];
416 Float_t q1 = clusterSignal[1];
417 Float_t q2 = clusterSignal[2];
418 Float_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
419 (clusterCharge*clusterCharge);
a819a5f7 420
421 // Correct for ExB displacement
17b26de4 422 if (fPar->ExBOn()) {
a819a5f7 423 Int_t local_time_bin = (Int_t) clusterPads[2];
424 Float_t driftLength = local_time_bin * timeBinSize + kAmWidth;
17b26de4 425 Float_t colSize = fPar->GetColPadSize(iplan);
a819a5f7 426 Float_t deltaY = omegaTau*driftLength/colSize;
427 clusterPads[1] = clusterPads[1] - deltaY;
428 }
429
47517f42 430 if (fVerbose > 1) {
3e1a3ad8 431 printf("-----------------------------------------------------------\n");
432 printf("Create cluster no. %d\n",nClusters);
433 printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0]
434 ,clusterPads[1]
435 ,clusterPads[2]);
436 printf("Indices: %d, %d, %d\n",clusterDigit[0]
437 ,clusterDigit[1]
438 ,clusterDigit[2]);
439 printf("Total charge = %f\n",clusterCharge);
440 printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0]
441 ,clusterTracks[1]
442 ,clusterTracks[2]);
443 printf(" pad1 %d, %d, %d\n",clusterTracks[3]
444 ,clusterTracks[4]
445 ,clusterTracks[5]);
446 printf(" pad2 %d, %d, %d\n",clusterTracks[6]
447 ,clusterTracks[7]
448 ,clusterTracks[8]);
db30bf0f 449 printf("Type = %d, Number of pads = %d\n",iType,nPadCount);
f7336fa3 450 }
451
5443e65e 452 // Calculate the position and the error
453 Float_t clusterPos[3];
454 clusterPos[0] = clusterPads[1] * colSize + col0;
455 clusterPos[1] = clusterPads[0] * rowSize + row0;
456 clusterPos[2] = clusterPads[2];
457 Float_t clusterSig[2];
458 clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
459 clusterSig[1] = rowSize * rowSize / 12.;
460
f7336fa3 461 // Add the cluster to the output array
bdbb05bb 462 AddCluster(clusterPos
463 ,idet
464 ,clusterCharge
465 ,clusterTracks
466 ,clusterSig
467 ,iType);
f7336fa3 468
469 }
3e1a3ad8 470 }
471 }
472 }
f7336fa3 473
3e1a3ad8 474 // Compress the arrays
475 digits->Compress(1,0);
476 track0->Compress(1,0);
477 track1->Compress(1,0);
478 track2->Compress(1,0);
f7336fa3 479
3e1a3ad8 480 // Write the cluster and reset the array
793ff80c 481 WriteClusters(idet);
bdbb05bb 482 ResetRecPoints();
793ff80c 483
47517f42 484 if (fVerbose > 0) {
17b26de4 485 printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42 486 printf("Found %d clusters in total.\n"
487 ,nClusters);
488 printf(" 2pad: %d\n",nClusters2pad);
489 printf(" 3pad: %d\n",nClusters3pad);
490 printf(" 4pad: %d\n",nClusters4pad);
491 printf(" 5pad: %d\n",nClusters5pad);
492 printf(" Large: %d\n",nClustersLarge);
493 }
f7336fa3 494
3e1a3ad8 495 }
496 }
497 }
f7336fa3 498
47517f42 499 if (fVerbose > 0) {
17b26de4 500 printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42 501 printf("Done.\n");
502 }
f7336fa3 503
504 return kTRUE;
505
506}
507
508//_____________________________________________________________________________
17b26de4 509Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Int_t plane, Float_t* padSignal)
f7336fa3 510{
511 //
512 // Method to unfold neighbouring maxima.
513 // The charge ratio on the overlapping pad is calculated
514 // until there is no more change within the range given by eps.
515 // The resulting ratio is then returned to the calling method.
516 //
517
17b26de4 518 Int_t irc = 0;
3e1a3ad8 519 Int_t itStep = 0; // Count iteration steps
f7336fa3 520
3e1a3ad8 521 Float_t ratio = 0.5; // Start value for ratio
522 Float_t prevRatio = 0; // Store previous ratio
f7336fa3 523
3e1a3ad8 524 Float_t newLeftSignal[3] = {0}; // Array to store left cluster signal
525 Float_t newRightSignal[3] = {0}; // Array to store right cluster signal
17b26de4 526 Float_t newSignal[3] = {0};
f7336fa3 527
3e1a3ad8 528 // Start the iteration
f7336fa3 529 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
530
531 itStep++;
532 prevRatio = ratio;
533
3e1a3ad8 534 // Cluster position according to charge ratio
535 Float_t maxLeft = (ratio*padSignal[2] - padSignal[0])
536 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
537 Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
538 / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
f7336fa3 539
3e1a3ad8 540 // Set cluster charge ratio
17b26de4 541 irc = fPar->PadResponse(1.0,maxLeft ,plane,newSignal);
542 Float_t ampLeft = padSignal[1] / newSignal[1];
543 irc = fPar->PadResponse(1.0,maxRight,plane,newSignal);
544 Float_t ampRight = padSignal[3] / newSignal[1];
f7336fa3 545
3e1a3ad8 546 // Apply pad response to parameters
17b26de4 547 irc = fPar->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
548 irc = fPar->PadResponse(ampRight,maxRight,plane,newRightSignal);
f7336fa3 549
3e1a3ad8 550 // Calculate new overlapping ratio
26edf6a4 551 ratio = TMath::Min((Float_t)1.0,newLeftSignal[2] /
db30bf0f 552 (newLeftSignal[2] + newRightSignal[0]));
f7336fa3 553
554 }
555
556 return ratio;
557
558}
559