ATO-123 - Configurin of the stream level. Swith to dump all signal withour zerro...
[u/mrichter/AliRoot.git] / TPC / TPCsim / AliTPCDigitizer.cxx
CommitLineData
3c038d07 1/**************************************************************************
2 * Copyright(c) 1998-2000, 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$ */
d37f8c34 17
f94f8f87 18/*
19 Class for creating of the sumable digits and digits from MC data
20 //
21 The input : ideal signals (Hits->Diffusion->Attachment -Ideal signal)
f9c05605 22 The output: "raw digits"
f94f8f87 23
24 Effect implemented:
25 1. Pad by pad gain map
f9c05605 26 2. Noise map
f94f8f87 27 3. The dead channels identified - zerro noise for corresponding pads
28 In this case the outpu equal zerro
f9c05605 29
f94f8f87 30*/
31
32
33
34
d37f8c34 35#include <stdlib.h>
3c038d07 36#include <TTree.h>
37#include <TObjArray.h>
38#include <TFile.h>
39#include <TDirectory.h>
19364939 40#include <Riostream.h>
f546a4b4 41#include <TParameter.h>
3c038d07 42
43#include "AliTPCDigitizer.h"
44
45#include "AliTPC.h"
46#include "AliTPCParam.h"
7a09f434 47#include "AliTPCParamSR.h"
3c038d07 48#include "AliRun.h"
c93255fe 49#include "AliLoader.h"
3c038d07 50#include "AliPDG.h"
f21fc003 51#include "AliDigitizationInput.h"
3c038d07 52#include "AliSimDigits.h"
f77f13c8 53#include "AliLog.h"
3c038d07 54
13116aec 55#include "AliTPCcalibDB.h"
56#include "AliTPCCalPad.h"
57#include "AliTPCCalROC.h"
0cbdc86b 58#include "TTreeStream.h"
f0254458 59#include "AliTPCReconstructor.h"
6fadf71b 60#include <TGraphErrors.h>
13116aec 61
a11596ad 62using std::cout;
63using std::cerr;
64using std::endl;
3c038d07 65ClassImp(AliTPCDigitizer)
66
67//___________________________________________
cd7d232c 68 AliTPCDigitizer::AliTPCDigitizer() :AliDigitizer(),fDebug(0), fDebugStreamer(0)
3c038d07 69{
179c6296 70 //
3c038d07 71// Default ctor - don't use it
179c6296 72//
73
3c038d07 74}
75
76//___________________________________________
f21fc003 77AliTPCDigitizer::AliTPCDigitizer(AliDigitizationInput* digInput)
cd7d232c 78 :AliDigitizer(digInput),fDebug(0), fDebugStreamer(0)
3c038d07 79{
179c6296 80 //
3c038d07 81// ctor which should be used
179c6296 82//
f21fc003 83 AliDebug(2,"(AliDigitizationInput* digInput) was processed");
d1570fe7 84 if (AliTPCReconstructor::StreamLevel()>0) fDebugStreamer = new TTreeSRedirector("TPCDigitDebug.root","recreate");
cd7d232c 85
3c038d07 86}
87
88//------------------------------------------------------------------------
89AliTPCDigitizer::~AliTPCDigitizer()
90{
91// Destructor
cd7d232c 92 if (fDebugStreamer) delete fDebugStreamer;
3c038d07 93}
94
95
96
97//------------------------------------------------------------------------
98Bool_t AliTPCDigitizer::Init()
99{
100// Initialization
101
102 return kTRUE;
103}
104
105
106//------------------------------------------------------------------------
f21fc003 107void AliTPCDigitizer::Digitize(Option_t* option)
3c038d07 108{
5b9e566f 109 //DigitizeFast(option);
110 DigitizeWithTailAndCrossTalk(option);
526c76a5 111
f648982e 112}
113//------------------------------------------------------------------------
f21fc003 114void AliTPCDigitizer::DigitizeFast(Option_t* option)
f648982e 115{
3c038d07 116
117 // merge input tree's with summable digits
118 //output stored in TreeTPCD
7a09f434 119 char s[100];
120 char ss[100];
3c038d07 121 TString optionString = option;
94bf739c 122 if (!strcmp(optionString.Data(),"deb")) {
f21fc003 123 cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
3c038d07 124 fDebug = 3;
125 }
126 //get detector and geometry
88cb7938 127
128
129 AliRunLoader *rl, *orl;
130 AliLoader *gime, *ogime;
131
132 if (gAlice == 0x0)
133 {
f21fc003 134 Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
135 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
88cb7938 136 if (rl == 0x0)
137 {
f21fc003 138 Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
88cb7938 139 return;
140 }
141 rl->LoadgAlice();
142 rl->GetAliRun();
143 }
3c038d07 144 AliTPC *pTPC = (AliTPC *) gAlice->GetModule("TPC");
145 AliTPCParam * param = pTPC->GetParam();
7a09f434 146
94e6c6f4 147 //sprintf(s,param->GetTitle());
f21fc003 148 snprintf(s,100,"%s",param->GetTitle());
94e6c6f4 149 //sprintf(ss,"75x40_100x60");
150 snprintf(ss,100,"75x40_100x60");
7a09f434 151 if(strcmp(s,ss)==0){
152 printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
153 delete param;
154 param=new AliTPCParamSR();
155 }
156 else{
94e6c6f4 157 //sprintf(ss,"75x40_100x60_150x60");
158 snprintf(ss,100,"75x40_100x60_150x60");
7a09f434 159 if(strcmp(s,ss)!=0) {
160 printf("No TPC parameters found...\n");
161 exit(2);
162 }
163 }
164
f546a4b4 165 pTPC->GenerNoise(500000); //create table with noise
407ff276 166 //
f21fc003 167 Int_t nInputs = fDigInput->GetNinputs();
407ff276 168 Int_t * masks = new Int_t[nInputs];
169 for (Int_t i=0; i<nInputs;i++)
f21fc003 170 masks[i]= fDigInput->GetMask(i);
407ff276 171 Short_t **pdig= new Short_t*[nInputs]; //pointers to the expanded digits array
f546a4b4 172 Int_t **ptr= new Int_t*[nInputs]; //pointers to the expanded tracks array
4df28b43 173 Bool_t *active= new Bool_t[nInputs]; //flag for active input segments
f546a4b4 174 Char_t phname[100];
88cb7938 175
3c038d07 176 //create digits array for given sectors
f648982e 177 // make indexes
f546a4b4 178 //
179 //create branch's in TPC treeD
f21fc003 180 orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
f546a4b4 181 ogime = orl->GetLoader("TPCLoader");
182 TTree * tree = ogime->TreeD();
183 AliSimDigits * digrow = new AliSimDigits;
184
185 if (tree == 0x0)
186 {
187 ogime->MakeTree("D");
188 tree = ogime->TreeD();
189 }
190 tree->Branch("Segment","AliSimDigits",&digrow);
191 //
407ff276 192 AliSimDigits ** digarr = new AliSimDigits*[nInputs];
88cb7938 193 for (Int_t i1=0;i1<nInputs; i1++)
194 {
195 digarr[i1]=0;
196 // intree[i1]
f21fc003 197 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
88cb7938 198 gime = rl->GetLoader("TPCLoader");
199 gime->LoadSDigits("read");
200 TTree * treear = gime->TreeS();
201
202 if (!treear)
203 {
204 cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
205 <<" input "<< i1<<endl;
712976a6 206 for (Int_t i2=0;i2<i1+1; i2++){
207
208 if(digarr[i2]) delete digarr[i2];
209 }
210 delete [] digarr;
211 delete [] active;
212 delete []masks;
213 delete []pdig;
214 delete []ptr;
88cb7938 215 return;
216 }
f546a4b4 217
94e6c6f4 218 //sprintf(phname,"lhcphase%d",i1);
219 snprintf(phname,100,"lhcphase%d",i1);
f546a4b4 220 TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
221 ->FindObject("lhcphase0");
222 if(!ph){
223 cerr<<"AliTPCDigitizer: LHC phase not found in"
224 <<" input "<< i1<<endl;
712976a6 225 for (Int_t i2=0;i2<i1+1; i2++){
226 if(digarr[i2]) delete digarr[i2];
227 }
228 delete [] digarr;
229 delete [] active;
230 delete []masks;
231 delete []pdig;
232 delete []ptr;
f546a4b4 233 return;
234 }
235 tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
236 //
88cb7938 237 if (treear->GetIndex()==0)
238 treear->BuildIndex("fSegmentID","fSegmentID");
239 treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
12d8d654 240 }
88cb7938 241
3c038d07 242
88cb7938 243
f546a4b4 244
3c038d07 245 //
246
407ff276 247 param->SetZeroSup(2);
248
f648982e 249 Int_t zerosup = param->GetZeroSup();
5033a39a 250 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
1ac191a6 251 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
f648982e 252 //
3c038d07 253 //Loop over segments of the TPC
254
4df28b43 255 for (Int_t segmentID=0; segmentID<param->GetNRowsTotal(); segmentID++)
88cb7938 256 {
f9c05605 257 Int_t sector, padRow;
258 if (!param->AdjustSectorRow(segmentID,sector,padRow))
4df28b43 259 {
260 cerr<<"AliTPC warning: invalid segment ID ! "<<segmentID<<endl;
261 continue;
262 }
f9c05605 263 AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector); // pad gains per given sector
264 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector); // noise per given sector
4df28b43 265 digrow->SetID(segmentID);
266
f9c05605 267 Int_t nTimeBins = 0;
268 Int_t nPads = 0;
4df28b43 269
270 Bool_t digitize = kFALSE;
271 for (Int_t i=0;i<nInputs; i++)
88cb7938 272 {
273
f21fc003 274 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
88cb7938 275 gime = rl->GetLoader("TPCLoader");
276
4df28b43 277 if (gime->TreeS()->GetEntryWithIndex(segmentID,segmentID) >= 0) {
278 digarr[i]->ExpandBuffer();
279 digarr[i]->ExpandTrackBuffer();
f9c05605 280 nTimeBins = digarr[i]->GetNRows();
281 nPads = digarr[i]->GetNCols();
4df28b43 282 active[i] = kTRUE;
f21fc003 283 if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
4df28b43 284 } else {
285 active[i] = kFALSE;
286 }
f21fc003 287 if (GetRegionOfInterest() && !digitize) break;
88cb7938 288 }
4df28b43 289 if (!digitize) continue;
3c038d07 290
f9c05605 291 digrow->Allocate(nTimeBins,nPads);
3c038d07 292 digrow->AllocateTrack(3);
293
407ff276 294 Float_t q=0;
295 Int_t label[1000]; //stack for 300 events
296 Int_t labptr = 0;
297
f9c05605 298 Int_t nElems = nTimeBins*nPads;
f648982e 299
88cb7938 300 for (Int_t i=0;i<nInputs; i++)
4df28b43 301 if (active[i]) {
88cb7938 302 pdig[i] = digarr[i]->GetDigits();
303 ptr[i] = digarr[i]->GetTracks();
4df28b43 304 }
88cb7938 305
407ff276 306 Short_t *pdig1= digrow->GetDigits();
307 Int_t *ptr1= digrow->GetTracks() ;
308
f648982e 309
407ff276 310
88cb7938 311 for (Int_t elem=0;elem<nElems; elem++)
312 {
313
314 q=0;
315 labptr=0;
316 // looop over digits
4df28b43 317 for (Int_t i=0;i<nInputs; i++) if (active[i])
88cb7938 318 {
319 // q += digarr[i]->GetDigitFast(rows,col);
320 q += *(pdig[i]);
321
322 for (Int_t tr=0;tr<3;tr++)
323 {
324 // Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
325 Int_t lab = ptr[i][tr*nElems];
326 if ( (lab > 1) && *(pdig[i])>zerosup)
327 {
328 label[labptr]=lab+masks[i];
329 labptr++;
330 }
331 }
332 pdig[i]++;
333 ptr[i]++;
334 }
335 q/=16.; //conversion factor
f9c05605 336 Float_t gain = gainROC->GetValue(padRow,elem/nTimeBins); // get gain for given - pad-row pad
ebf9a79e 337 //if (gain<0.5){
338 //printf("problem\n");
339 //}
13116aec 340 q*= gain;
f9c05605 341 Float_t noisePad = noiseROC->GetValue(padRow,elem/nTimeBins);
88cb7938 342 // Float_t noise = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
343 Float_t noise = pTPC->GetNoise();
1ac191a6 344 q+=noise*noisePad;
3c038d07 345 q=TMath::Nint(q);
88cb7938 346 if (q > zerosup)
347 {
e93a083a 348 if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
88cb7938 349 //digrow->SetDigitFast((Short_t)q,rows,col);
350 *pdig1 =Short_t(q);
351 for (Int_t tr=0;tr<3;tr++)
352 {
353 if (tr<labptr)
354 ptr1[tr*nElems] = label[tr];
355 }
356 }
357 pdig1++;
358 ptr1++;
359 }
58c45a5e 360 //
361 // glitch filter
362 //
363 digrow->GlitchFilter();
364 //
3c038d07 365 digrow->CompresBuffer(1,zerosup);
366 digrow->CompresTrackBuffer(1);
367 tree->Fill();
f9c05605 368 if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
4df28b43 369 } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
88cb7938 370
371
f21fc003 372 orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
88cb7938 373 ogime = orl->GetLoader("TPCLoader");
374 ogime->WriteDigits("OVERWRITE");
375
f21fc003 376 //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
88cb7938 377
3c038d07 378 delete digrow;
407ff276 379 for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
380 delete []masks;
4df28b43 381 delete []pdig;
382 delete []ptr;
383 delete []active;
384 delete []digarr;
3c038d07 385}
386
f648982e 387
388
389//------------------------------------------------------------------------
f21fc003 390void AliTPCDigitizer::DigitizeSave(Option_t* option)
f648982e 391{
f9c05605 392 //
393 // Merge input tree's with summable digits
394 // Output digits stored in TreeTPCD
395 //
396 // Not active for long time.
397 // Before adding modification (for ion tail calucation and for the crorsstalk) it should be
398 // checked one by one with currenlty used AliTPCDigitizer::DigitizeFast
399 //
f648982e 400 TString optionString = option;
94bf739c 401 if (!strcmp(optionString.Data(),"deb")) {
f21fc003 402 cout<<"AliTPCDigitizer::Digitize: called with option deb "<<endl;
f648982e 403 fDebug = 3;
404 }
405 //get detector and geometry
88cb7938 406 AliRunLoader *rl, *orl;
407 AliLoader *gime, *ogime;
408
409
f21fc003 410 orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
88cb7938 411 ogime = orl->GetLoader("TPCLoader");
412
f21fc003 413 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
097d753a 414 //gime = rl->GetLoader("TPCLoader");
415 rl->GetLoader("TPCLoader");
88cb7938 416 rl->LoadgAlice();
417 AliRun* alirun = rl->GetAliRun();
418
419 AliTPC *pTPC = (AliTPC *) alirun->GetModule("TPC");
f648982e 420 AliTPCParam * param = pTPC->GetParam();
421 pTPC->GenerNoise(500000); //create teble with noise
422 printf("noise %f \n", param->GetNoise()*param->GetNoiseNormFac());
423 //
f21fc003 424 Int_t nInputs = fDigInput->GetNinputs();
097d753a 425 // stupid protection...
426 if (nInputs <= 0) return;
427 //
f648982e 428 Int_t * masks = new Int_t[nInputs];
429 for (Int_t i=0; i<nInputs;i++)
f21fc003 430 masks[i]= fDigInput->GetMask(i);
f648982e 431
606d6545 432 AliSimDigits ** digarr = new AliSimDigits*[nInputs];
433 for(Int_t ii=0;ii<nInputs;ii++) digarr[ii]=0;
434
88cb7938 435 for (Int_t i1=0;i1<nInputs; i1++)
436 {
606d6545 437 //digarr[i1]=0;
f648982e 438 // intree[i1]
f21fc003 439 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
88cb7938 440 gime = rl->GetLoader("TPCLoader");
441
442 TTree * treear = gime->TreeS();
de0e19ab 443 //
f648982e 444 if (!treear) {
de0e19ab 445 cerr<<" TPC - not existing input = \n"<<i1<<" ";
446 delete [] masks;
447 for(Int_t i=0; i<nInputs; i++) delete digarr[i];
448 delete [] digarr;
449 return;
f648982e 450 }
de0e19ab 451 //
452 TBranch * br = treear->GetBranch("fSegmentID");
453 if (br) br->GetFile()->cd();
f648982e 454 treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
de0e19ab 455 }
88cb7938 456
f21fc003 457 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
88cb7938 458 gime = rl->GetLoader("TPCLoader");
459 Stat_t nentries = gime->TreeS()->GetEntries();
f648982e 460
461
462 //create branch's in TPC treeD
463 AliSimDigits * digrow = new AliSimDigits;
88cb7938 464 TTree * tree = ogime->TreeD();
465
f648982e 466 tree->Branch("Segment","AliSimDigits",&digrow);
f648982e 467 param->SetZeroSup(2);
468
469 Int_t zerosup = param->GetZeroSup();
470 //Loop over segments of the TPC
471
5033a39a 472 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
1ac191a6 473 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
f648982e 474 for (Int_t n=0; n<nentries; n++) {
f21fc003 475 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
88cb7938 476 gime = rl->GetLoader("TPCLoader");
477 gime->TreeS()->GetEvent(n);
478
f648982e 479 digarr[0]->ExpandBuffer();
480 digarr[0]->ExpandTrackBuffer();
13116aec 481
482
f648982e 483 for (Int_t i=1;i<nInputs; i++){
f21fc003 484// fDigInput->GetInputTreeTPCS(i)->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
485 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
88cb7938 486 gime = rl->GetLoader("TPCLoader");
487 gime->TreeS()->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
f648982e 488 digarr[i]->ExpandBuffer();
489 digarr[i]->ExpandTrackBuffer();
490 if ((digarr[0]->GetID()-digarr[i]->GetID())>0)
88cb7938 491 printf("problem\n");
f648982e 492
493 }
494
f9c05605 495 Int_t sector, padRow;
496 if (!param->AdjustSectorRow(digarr[0]->GetID(),sector,padRow)) {
f648982e 497 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
498 continue;
499 }
500
f9c05605 501 AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector); // pad gains per given sector
502 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector); // noise per given sector
f648982e 503 digrow->SetID(digarr[0]->GetID());
504
f9c05605 505 Int_t nTimeBins = digarr[0]->GetNRows();
506 Int_t nPads = digarr[0]->GetNCols();
507 digrow->Allocate(nTimeBins,nPads);
f648982e 508 digrow->AllocateTrack(3);
509
510 Float_t q=0;
511 Int_t label[1000]; //stack for 300 events
512 Int_t labptr = 0;
513
514
515
f9c05605 516 for (Int_t iTimeBin=0;iTimeBin<nTimeBins; iTimeBin++){ // iTimeBin
517 for (Int_t iPad=0;iPad<nPads; iPad++){ // pad
f648982e 518
88cb7938 519 q=0;
520 labptr=0;
521 // looop over digits
f648982e 522 for (Int_t i=0;i<nInputs; i++){
f9c05605 523 q += digarr[i]->GetDigitFast(iTimeBin,iPad);
f648982e 524 //q += *(pdig[i]);
88cb7938 525
f648982e 526 for (Int_t tr=0;tr<3;tr++) {
f9c05605 527 Int_t lab = digarr[i]->GetTrackIDFast(iTimeBin,iPad,tr);
88cb7938 528 //Int_t lab = ptr[i][tr*nElems];
f648982e 529 if ( (lab > 1) ) {
530 label[labptr]=lab+masks[i];
531 labptr++;
88cb7938 532 }
f648982e 533 }
88cb7938 534 // pdig[i]++;
535 //ptr[i]++;
536
f648982e 537 }
88cb7938 538 q/=16.; //conversion factor
539 // Float_t noise = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
f9c05605 540 Float_t gain = gainROC->GetValue(padRow,iPad);
13116aec 541 q*= gain;
f9c05605 542 Float_t noisePad = noiseROC->GetValue(padRow, iPad);
1ac191a6 543
88cb7938 544 Float_t noise = pTPC->GetNoise();
1ac191a6 545 q+=noise*noisePad;
f9c05605 546 //
547 // here we can get digits from past and add signal
548 //
549 //
550 //for (Int_t jTimeBin=0; jTimeBin<iTimeBin; jTimeBin++)
551 // q+=ionTail
552 //
633fa37f 553
f648982e 554 q=TMath::Nint(q);
555 if (q > zerosup){
88cb7938 556
e93a083a 557 if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
f9c05605 558 digrow->SetDigitFast((Short_t)q,iTimeBin,iPad);
88cb7938 559 // *pdig1 =Short_t(q);
560 for (Int_t tr=0;tr<3;tr++){
561 if (tr<labptr)
f9c05605 562 ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],iTimeBin,iPad,tr);
88cb7938 563 //ptr1[tr*nElems] = label[tr];
564 //else
f9c05605 565 // ((AliSimDigits*)digrow)->SetTrackIDFast(-1,iTimeBin,iPad,tr);
88cb7938 566 // ptr1[tr*nElems] = 1;
567 }
568 }
569 //pdig1++;
570 //ptr1++;
f648982e 571 }
572 }
573
574 digrow->CompresBuffer(1,zerosup);
575 digrow->CompresTrackBuffer(1);
576 tree->Fill();
f9c05605 577 if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
f648982e 578 }
f21fc003 579// printf("end TPC merging - end -Tree %s\t%p\n",fDigInput->GetInputTreeH(0)->GetName(),fDigInput->GetInputTreeH(0)->GetListOfBranches()->At(3));
580 //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
88cb7938 581 ogime->WriteDigits("OVERWRITE");
582
583 for (Int_t i=1;i<nInputs; i++)
584 {
f21fc003 585 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
88cb7938 586 gime = rl->GetLoader("TPCLoader");
587 gime->UnloadSDigits();
588 }
589 ogime->UnloadDigits();
590
f648982e 591 delete digrow;
592 for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
4ce766eb 593 delete [] masks;
594 delete [] digarr;
f648982e 595}
0cbdc86b 596
597
598
6fadf71b 599
600
601
602
0cbdc86b 603//------------------------------------------------------------------------
cd7d232c 604void AliTPCDigitizer::DigitizeWithTailAndCrossTalk(Option_t* option)
0cbdc86b 605{
606 // Modified version of the digitization function
607 // Modification: adding the ion tail and crosstalk:
608 //
609 // pcstream used in order to visually inspect data
610 //
611 //
612 // Crosstalk simulation:
613 // 1.) Calculate per time bin mean charge (per pad) within anode wire segment
614 // 2.) Subsract for the clusters at given time bin fraction of (mean charge) normalized by add hoc constant
615 // AliTPCRecoParam::GetCrosstalkCorrection() (0 if not crosstalk, 1 if ideal crosstalk)
616 // for simplicity we are assuming that wire segents are related to pad-rows
617 // Wire segmentationn is obtatined from the
618 // AliTPCParam::GetWireSegment(Int_t sector, Int_t row); // to be implemented
619 // AliTPCParam::GetNPadsPerSegment(Int_t segmentID); // to be implemented
f5a7275d 620 // 3.) Substract form the signal contribution from the previous tracks - Ion tail in case speified in the AliTPCRecoParam
621 // AliTPCRecoParam::GetUseIonTailCorrection()
0cbdc86b 622 //
623 // Ion tail simulation:
624 // 1.) Needs signal from pad+-1, taking signal from history
625 // merge input tree's with summable digits
626 // output stored in TreeTPCD
f5a7275d 627 //
628 AliTPCcalibDB* const calib=AliTPCcalibDB::Instance();
629 AliTPCRecoParam *recoParam = calib->GetRecoParam(0);
d1570fe7 630 AliDebug(1, Form(" recoParam->GetCrosstalkCorrection() = %f", recoParam->GetCrosstalkCorrection()));
dbc380a1 631 AliDebug(1,Form(" recoParam->GetUseIonTailCorrection() = %d ", recoParam->GetUseIonTailCorrection()));
6fadf71b 632 Int_t nROCs = 72;
0cbdc86b 633 char s[100];
634 char ss[100];
635 TString optionString = option;
636 if (!strcmp(optionString.Data(),"deb")) {
637 cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
638 fDebug = 3;
639 }
0cbdc86b 640
6fadf71b 641 // ======== get detector and geometry =======
0cbdc86b 642
643 AliRunLoader *rl, *orl;
644 AliLoader *gime, *ogime;
6fadf71b 645
0cbdc86b 646 if (gAlice == 0x0)
6fadf71b 647 {
648 Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
649 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
650 if (rl == 0x0)
651 {
652 Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
653 return;
654 }
655 rl->LoadgAlice();
656 rl->GetAliRun();
657 }
0cbdc86b 658 AliTPC *pTPC = (AliTPC *) gAlice->GetModule("TPC");
659 AliTPCParam * param = pTPC->GetParam();
6fadf71b 660
0cbdc86b 661 //sprintf(s,param->GetTitle());
662 snprintf(s,100,"%s",param->GetTitle());
663 //sprintf(ss,"75x40_100x60");
664 snprintf(ss,100,"75x40_100x60");
665 if(strcmp(s,ss)==0){
666 printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
667 delete param;
668 param=new AliTPCParamSR();
6fadf71b 669 } else {
0cbdc86b 670 //sprintf(ss,"75x40_100x60_150x60");
671 snprintf(ss,100,"75x40_100x60_150x60");
6fadf71b 672 if(strcmp(s,ss)!=0) {
673 printf("No TPC parameters found...\n");
674 exit(2);
675 }
0cbdc86b 676 }
6fadf71b 677
0cbdc86b 678 pTPC->GenerNoise(500000); //create table with noise
679 //
680 Int_t nInputs = fDigInput->GetNinputs();
681 Int_t * masks = new Int_t[nInputs];
682 for (Int_t i=0; i<nInputs;i++)
683 masks[i]= fDigInput->GetMask(i);
684 Short_t **pdig= new Short_t*[nInputs]; //pointers to the expanded digits array
685 Int_t **ptr= new Int_t*[nInputs]; //pointers to the expanded tracks array
686 Bool_t *active= new Bool_t[nInputs]; //flag for active input segments
687 Char_t phname[100];
6fadf71b 688
0cbdc86b 689 //create digits array for given sectors
690 // make indexes
691 //
692 //create branch's in TPC treeD
693 orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
694 ogime = orl->GetLoader("TPCLoader");
695 TTree * tree = ogime->TreeD();
696 AliSimDigits * digrow = new AliSimDigits;
697
698 if (tree == 0x0)
6fadf71b 699 {
700 ogime->MakeTree("D");
701 tree = ogime->TreeD();
702 }
0cbdc86b 703 tree->Branch("Segment","AliSimDigits",&digrow);
704 //
705 AliSimDigits ** digarr = new AliSimDigits*[nInputs];
706 for (Int_t i1=0;i1<nInputs; i1++)
6fadf71b 707 {
708 digarr[i1]=0;
709 // intree[i1]
710 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
711 gime = rl->GetLoader("TPCLoader");
712 gime->LoadSDigits("read");
713 TTree * treear = gime->TreeS();
714
715 if (!treear)
0cbdc86b 716 {
6fadf71b 717 cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
718 <<" input "<< i1<<endl;
719 for (Int_t i2=0;i2<i1+1; i2++){
0cbdc86b 720
6fadf71b 721 if(digarr[i2]) delete digarr[i2];
0cbdc86b 722 }
6fadf71b 723 delete [] digarr;
724 delete [] active;
725 delete []masks;
726 delete []pdig;
727 delete []ptr;
728 return;
0cbdc86b 729 }
730
6fadf71b 731 //sprintf(phname,"lhcphase%d",i1);
732 snprintf(phname,100,"lhcphase%d",i1);
733 TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
734 ->FindObject("lhcphase0");
735 if(!ph)
736 {
737 cerr<<"AliTPCDigitizer: LHC phase not found in"
738 <<" input "<< i1<<endl;
739 for (Int_t i2=0;i2<i1+1; i2++){
740 if(digarr[i2]) delete digarr[i2];
741 }
742 delete [] digarr;
743 delete [] active;
744 delete []masks;
745 delete []pdig;
746 delete []ptr;
747 return;
748 }
749 tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
750 //
751 if (treear->GetIndex()==0)
752 treear->BuildIndex("fSegmentID","fSegmentID");
753 treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
754 }
0cbdc86b 755
756
757
0cbdc86b 758
6fadf71b 759 //
760 // zero supp, take gain and noise map of TPC from OCDB
0cbdc86b 761 param->SetZeroSup(2);
762 Int_t zerosup = param->GetZeroSup();
763 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
764 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
6fadf71b 765
766
767 //
768 // Cache the ion tail objects form OCDB
769 //
770
526c76a5 771 TObjArray *ionTailArr = (TObjArray*)AliTPCcalibDB::Instance()->GetIonTailArray();
772 if (!ionTailArr) {AliFatal("TPC - Missing IonTail OCDB object");}
0f155c90 773 // TObject *rocFactorIROC = ionTailArr->FindObject("factorIROC");
774// TObject *rocFactorOROC = ionTailArr->FindObject("factorOROC");
775// Float_t factorIROC = (atof(rocFactorIROC->GetTitle()));
776// Float_t factorOROC = (atof(rocFactorOROC->GetTitle()));
8cc1870a 777 Int_t nIonTailBins =0;
6fadf71b 778 TObjArray timeResFunc(nROCs);
779 for (Int_t isec = 0;isec<nROCs;isec++){ //loop overs sectors
6fadf71b 780 // Array of TGraphErrors for a given sector
781 TGraphErrors ** graphRes = new TGraphErrors *[20];
a848abbf 782 Float_t * trfIndexArr = new Float_t[20];
6fadf71b 783 for (Int_t icache=0; icache<20; icache++)
784 {
785 graphRes[icache] = NULL;
a848abbf 786 trfIndexArr[icache] = 0;
6fadf71b 787 }
a848abbf 788 if (!AliTPCcalibDB::Instance()->GetTailcancelationGraphs(isec,graphRes,trfIndexArr)) continue;
789
526c76a5 790 // fill all TGraphErrors of trfs (time response function) of a given sector to a TObjArray
6fadf71b 791 TObjArray *timeResArr = new TObjArray(20); timeResArr -> SetOwner(kTRUE);
792 for (Int_t ires = 0;ires<20;ires++) timeResArr->AddAt(graphRes[ires],ires);
526c76a5 793 timeResFunc.AddAt(timeResArr,isec); // Fill all trfs into a single TObjArray
8cc1870a 794 nIonTailBins = graphRes[3]->GetN();
6fadf71b 795 }
796
0cbdc86b 797 //
6fadf71b 798 // 1.) Make first loop to calculate mean amplitude per pad per segment for cross talk
0cbdc86b 799 //
6fadf71b 800
801 TObjArray crossTalkSignalArray(nROCs); // for 36 sectors
802 TVectorD * qTotSectorOld = new TVectorD(nROCs);
6fadf71b 803 Float_t qTotTPC = 0.;
d734bfaa 804 Float_t qTotPerSector = 0.;
6fadf71b 805 Int_t nTimeBinsAll = 1100;
806 Int_t nWireSegments=11;
807 // 1.a) crorstalk matrix initialization
808 for (Int_t sector=0; sector<nROCs; sector++){
809 TMatrixD *pcrossTalkSignal = new TMatrixD(nWireSegments,nTimeBinsAll);
810 for (Int_t imatrix = 0; imatrix<11; imatrix++)
811 for (Int_t jmatrix = 0; jmatrix<nTimeBinsAll; jmatrix++){
812 (*pcrossTalkSignal)[imatrix][jmatrix]=0.;
813 }
814 crossTalkSignalArray.AddAt(pcrossTalkSignal,sector);
815 }
816 //
817 // main loop over rows of whole TPC
818 for (Int_t globalRowID=0; globalRowID<param->GetNRowsTotal(); globalRowID++) {
819 Int_t sector, padRow;
820 if (!param->AdjustSectorRow(globalRowID,sector,padRow)) {
821 cerr<<"AliTPC warning: invalid segment ID ! "<<globalRowID<<endl;
822 continue;
823 }
824 // Calculate number of pads in a anode wire segment for normalization
d734bfaa 825 Int_t wireSegmentID = param->GetWireSegment(sector,padRow);
826 Float_t nPadsPerSegment = (Float_t)(param->GetNPadsPerSegment(wireSegmentID));
6fadf71b 827 // structure with mean signal per pad to be filled for each timebin in first loop (11 anodeWireSegment and 1100 timebin)
828 TMatrixD &crossTalkSignal = *((TMatrixD*)crossTalkSignalArray.At(sector));
829 AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector); // pad gains per given sector
d1570fe7 830 // digrow->SetID(globalRowID);
6fadf71b 831 Int_t nTimeBins = 0;
832 Int_t nPads = 0;
833 Bool_t digitize = kFALSE;
834 for (Int_t i=0;i<nInputs; i++){ //here we can have more than one input - merging of separate events, signal1+signal2+background
835 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
836 gime = rl->GetLoader("TPCLoader");
837 if (gime->TreeS()->GetEntryWithIndex(globalRowID,globalRowID) >= 0) {
838 digarr[i]->ExpandBuffer();
839 digarr[i]->ExpandTrackBuffer();
840 nTimeBins = digarr[i]->GetNRows();
841 nPads = digarr[i]->GetNCols();
842 active[i] = kTRUE;
843 if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
844 } else {
845 active[i] = kFALSE;
846 }
847 if (GetRegionOfInterest() && !digitize) break;
848 }
849 if (!digitize) continue;
d1570fe7 850 //digrow->Allocate(nTimeBins,nPads);
526c76a5 851 Float_t q = 0;
6fadf71b 852 Int_t labptr = 0;
853 Int_t nElems = nTimeBins*nPads; // element is a unit of a given row's pad-timebin space
854 for (Int_t i=0;i<nInputs; i++)
855 if (active[i]) {
856 pdig[i] = digarr[i]->GetDigits();
526c76a5 857 }
858 //
859 // loop over elements i.e pad-timebin space of a row, "padNumber=elem/nTimeBins", "timeBin=elem%nTimeBins"
860 //
6fadf71b 861 for (Int_t elem=0;elem<nElems; elem++) {
862 q=0;
863 labptr=0;
864 // looop over digits
865 for (Int_t i=0;i<nInputs; i++) if (active[i]) {
866 q += *(pdig[i]);
867 pdig[i]++;
868 }
d1570fe7 869 if (q<=0) continue;
6fadf71b 870 Int_t padNumber = elem/nTimeBins;
871 Int_t timeBin = elem%nTimeBins;
6fadf71b 872 Float_t gain = gainROC->GetValue(padRow,padNumber); // get gain for given - pad-row pad
873 q*= gain;
d734bfaa 874 crossTalkSignal[wireSegmentID][timeBin]+= q/nPadsPerSegment; // Qtot per segment for a given timebin
526c76a5 875 qTotSectorOld -> GetMatrixArray()[sector] += q; // Qtot for each sector
876 qTotTPC += q; // Qtot for whole TPC
6fadf71b 877 } // end of q loop
878 } // end of global row loop
879
0f155c90 880 //
881 // 1.b) Dump the content of the crossTalk signal to the debug stremer - to be corealted later with the same crosstalk correction
882 // assumed during reconstruction
883 //
3ec92887 884 if ((AliTPCReconstructor::StreamLevel()&kStreamCrosstalk)>0) {
0f155c90 885 //
886 // dump the crosstalk matrices to tree for further investigation
887 // a.) to estimate fluctuation of pedestal in indiviula wire segments
888 // b.) to check correlation between regions
889 // c.) to check relative conribution of signal below threshold to crosstalk
890 for (Int_t isector=0; isector<nROCs; isector++){ //set all ellemts of crosstalk matrix to 0
891 TMatrixD * crossTalkMatrix = (TMatrixD*)crossTalkSignalArray.At(isector);
892 //TMatrixD * crossTalkMatrixBelow = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs);
893 TVectorD vecAll(crossTalkMatrix->GetNrows());
894 //TVectorD vecBelow(crossTalkMatrix->GetNrows());
895 //
896 for (Int_t itime=0; itime<crossTalkMatrix->GetNcols(); itime++){
897 for (Int_t iwire=0; iwire<crossTalkMatrix->GetNrows(); iwire++){
898 vecAll[iwire]=(*crossTalkMatrix)(iwire,itime);
899 //vecBelow[iwire]=(*crossTalkMatrixBelow)(iwire,itime);
900 }
901 (*fDebugStreamer)<<"crosstalkMatrix"<<
902 "sector="<<isector<<
903 "itime="<<itime<<
904 "vecAll.="<<&vecAll<< // crosstalk charge + charge below threshold
905 //"vecBelow.="<<&vecBelow<< // crosstalk contribution from signal below threshold
906 "\n";
907 }
908 }
909 }
910
911
912
6fadf71b 913
0cbdc86b 914 //
6fadf71b 915 // 2.) Loop over segments (padrows) of the TPC
0cbdc86b 916 //
6fadf71b 917 for (Int_t globalRowID=0; globalRowID<param->GetNRowsTotal(); globalRowID++) {
0cbdc86b 918 Int_t sector, padRow;
6fadf71b 919 if (!param->AdjustSectorRow(globalRowID,sector,padRow)) {
920 cerr<<"AliTPC warning: invalid segment ID ! "<<globalRowID<<endl;
0cbdc86b 921 continue;
6fadf71b 922 }
a848abbf 923 TObjArray *arrTRF = (TObjArray*)timeResFunc.At(sector);
faec07d1 924 // TGraphErrors *graphTRF = (TGraphErrors*)arrTRF->At(1);
d734bfaa 925 Int_t wireSegmentID = param->GetWireSegment(sector,padRow);
926 Float_t nPadsPerSegment = (Float_t)(param->GetNPadsPerSegment(wireSegmentID));
faec07d1 927 // const Float_t ampfactor = (sector<36)?factorIROC:factorOROC; // factor for the iontail which is ROC type dependent
526c76a5 928 AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector); // pad gains per given sector
0cbdc86b 929 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector); // noise per given sector
6fadf71b 930 digrow->SetID(globalRowID);
0cbdc86b 931 Int_t nTimeBins = 0;
932 Int_t nPads = 0;
0cbdc86b 933 Bool_t digitize = kFALSE;
526c76a5 934 for (Int_t i=0;i<nInputs; i++) { //here we can have more than one input - merging of separate events, signal1+signal2+background
0cbdc86b 935 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
936 gime = rl->GetLoader("TPCLoader");
6fadf71b 937 if (gime->TreeS()->GetEntryWithIndex(globalRowID,globalRowID) >= 0) {
938 digarr[i]->ExpandBuffer();
939 digarr[i]->ExpandTrackBuffer();
0cbdc86b 940 nTimeBins = digarr[i]->GetNRows();
941 nPads = digarr[i]->GetNCols();
6fadf71b 942 active[i] = kTRUE;
943 if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
0cbdc86b 944 } else {
6fadf71b 945 active[i] = kFALSE;
0cbdc86b 946 }
947 if (GetRegionOfInterest() && !digitize) break;
6fadf71b 948 }
0cbdc86b 949 if (!digitize) continue;
6fadf71b 950
0cbdc86b 951 digrow->Allocate(nTimeBins,nPads);
952 digrow->AllocateTrack(3);
6fadf71b 953
d734bfaa 954 Int_t localPad = 0;
526c76a5 955 Float_t q = 0.;
d734bfaa 956 Float_t qXtalk = 0.;
957 Float_t qIonTail = 0.;
958 Float_t qOrig = 0.;
0cbdc86b 959 Int_t label[1000]; //stack for 300 events
960 Int_t labptr = 0;
6fadf71b 961 Int_t nElems = nTimeBins*nPads; // element is a unit of a given row's pad-timebin space
0cbdc86b 962 for (Int_t i=0;i<nInputs; i++)
6fadf71b 963 if (active[i]) {
964 pdig[i] = digarr[i]->GetDigits();
965 ptr[i] = digarr[i]->GetTracks();
0cbdc86b 966 }
0cbdc86b 967 Short_t *pdig1= digrow->GetDigits();
968 Int_t *ptr1= digrow->GetTracks() ;
6fadf71b 969 // loop over elements i.e pad-timebin space of a row
1c29da5d 970 for (Int_t elem=0;elem<nElems; elem++) {
526c76a5 971 q=0;
6fadf71b 972 labptr=0;
973 // looop over digits
1c29da5d 974 for (Int_t i=0;i<nInputs; i++) if (active[i]){
6fadf71b 975 q += *(pdig[i]);
1c29da5d 976 for (Int_t tr=0;tr<3;tr++) {
6fadf71b 977 Int_t lab = ptr[i][tr*nElems];
1c29da5d 978 if ( (lab > 1) && *(pdig[i])>zerosup) {
6fadf71b 979 label[labptr]=lab+masks[i];
980 labptr++;
981 }
982 }
983 pdig[i]++;
984 ptr[i]++;
985 }
526c76a5 986 Int_t padNumber = elem/nTimeBins;
987 Int_t timeBin = elem%nTimeBins;
d734bfaa 988 localPad = padNumber-nPads/2;
1c29da5d 989
6fadf71b 990 Float_t gain = gainROC->GetValue(padRow,padNumber); // get gain for given - pad-row pad
991 //if (gain<0.5){
992 //printf("problem\n");
993 //}
994 q*= gain;
d734bfaa 995 qOrig = q;
d1570fe7 996 Float_t noisePad = noiseROC->GetValue(padRow,padNumber);
997 Float_t noise = pTPC->GetNoise()*noisePad;
c4dc8a5b 998 if ( (q/16.+noise)> zerosup){
999 // Crosstalk correction
1000 qXtalk = (*(TMatrixD*)crossTalkSignalArray.At(sector))[wireSegmentID][timeBin];
1001 qTotPerSector = qTotSectorOld -> GetMatrixArray()[sector];
1002
1003 // Ion tail correction: being elem=padNumber*nTimeBins+timeBin;
1004 Int_t lowerElem=elem-nIonTailBins;
1005 Int_t zeroElem =(elem/nTimeBins)*nTimeBins;
1006 if (zeroElem<0) zeroElem=0;
1007 if (lowerElem<zeroElem) lowerElem=zeroElem;
1008 //
1009 qIonTail=0;
1010 if (q>0 && recoParam->GetUseIonTailCorrection()){
dbc380a1 1011 for (Int_t i=0;i<nInputs; i++) if (active[i]){
c4dc8a5b 1012 Short_t *pdigC= digarr[i]->GetDigits();
faec07d1 1013 if (padNumber==0) continue;
1014 if (padNumber>=nPads-1) continue;
1015 for (Int_t dpad=-1; dpad<=1; dpad++){ // calculate iontail due signals from neigborhood pads
dbc380a1 1016 for (Int_t celem=elem-1; celem>lowerElem; celem--){
1017 Int_t celemPad=celem+dpad*nTimeBins;
faec07d1 1018 Double_t qCElem=pdigC[celemPad];
1019 if ( qCElem<=0) continue;
dbc380a1 1020 //here we substract ion tail
1021 Double_t COG=0;
1022 if (celemPad-nTimeBins>nTimeBins && celemPad+nTimeBins<nElems){ // COG calculation in respect to current pad in pad units
1023 Double_t sumAmp=pdigC[celemPad-nTimeBins]+pdigC[celemPad]+pdigC[celemPad+nTimeBins];
1024 COG=(-1.0*pdigC[celemPad-nTimeBins]+pdigC[celemPad+nTimeBins])/sumAmp;
1025 }
faec07d1 1026 Int_t indexTRFPRF = (TMath::Nint(TMath::Abs(COG*10.))%20);
1027 TGraphErrors *graphTRFPRF = (TGraphErrors*)arrTRF->At(indexTRFPRF);
1028 if (graphTRFPRF==NULL) continue;
dbc380a1 1029 // here we should get index and point of TRF corresponding to given COG position
faec07d1 1030 if (graphTRFPRF->GetY()[elem-celem]<0)qIonTail+=qCElem*graphTRFPRF->GetY()[elem-celem];
dbc380a1 1031 }
1c29da5d 1032 }
c4dc8a5b 1033 }
1c29da5d 1034 }
1035 }
6fadf71b 1036 //
d1570fe7 1037 q -= qXtalk*recoParam->GetCrosstalkCorrection();
a848abbf 1038 q+=qIonTail;
1c29da5d 1039 q/=16.; //conversion factor
d1570fe7 1040 q+=noise;
6fadf71b 1041 q=TMath::Nint(q); // round to the nearest integer
526c76a5 1042
1043
1044 // fill info for degugging
3ec92887 1045 if ( ((AliTPCReconstructor::StreamLevel()&kStreamSignal)>0) && ((qOrig > zerosup)||((AliTPCReconstructor::StreamLevel()&kStreamSignalAll)>0) )) {
1046 TTreeSRedirector &cstream = *fDebugStreamer;
1047 UInt_t uid = AliTPCROC::GetTPCUniqueID(sector, padRow, padNumber);
526c76a5 1048 cstream <<"ionTailXtalk"<<
f618a956 1049 "uid="<<uid<< // globla unique identifier
d734bfaa 1050 "sector="<< sector<<
1051 "globalRowID="<<globalRowID<<
1052 "padRow="<< padRow<< //pad row
1053 "wireSegmentID="<< wireSegmentID<< //wire segment 0-11, 0-3 in IROC 4-10 in OROC
1054 "localPad="<<localPad<< // pad number -npads/2 .. npads/2
1055 "padNumber="<<padNumber<< // pad number 0..npads
1056 "timeBin="<< timeBin<< // time bin
1057 "nPadsPerSegment="<<nPadsPerSegment<<// number of pads per wire segment
1058 "qTotPerSector="<<qTotPerSector<< // total charge in sector
1059 //
1060 "qTotTPC="<<qTotTPC<< // acumulated charge without crosstalk and ion tail in full TPC
1061 "qOrig="<< qOrig<< // charge in given pad-row,pad,time-bin
1062 "q="<<q<< // q=qOrig-qXtalk-qIonTail - to check sign of the effects
1063 "qXtalk="<<qXtalk<< // crosstal contribtion at given position
1064 "qIonTail="<<qIonTail<< // ion tail cotribution from past signal
526c76a5 1065 "\n";
1066 }// dump the results to the debug streamer if in debug mode
1067
1068 if (q > zerosup){
6fadf71b 1069 if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
1070 //digrow->SetDigitFast((Short_t)q,rows,col);
1071 *pdig1 =Short_t(q);
1072 for (Int_t tr=0;tr<3;tr++)
526c76a5 1073 {
6fadf71b 1074 if (tr<labptr)
1075 ptr1[tr*nElems] = label[tr];
1076 }
1077 }
1078 pdig1++;
1079 ptr1++;
1080 }
d734bfaa 1081
8cc1870a 1082 if (AliTPCReconstructor::StreamLevel()==1 && qOrig > zerosup ) {
d734bfaa 1083 cout << " sector = " << sector << " row = " << padRow << " localPad = " << localPad << " SegmentID = " << wireSegmentID << " nPadsPerSegment = " << nPadsPerSegment << endl;
1084 cout << " qXtalk = " << qXtalk ;
1085 cout << " qOrig = " << qOrig ;
1086 cout << " q = " << q ;
1087 cout << " qsec = " << qTotPerSector ;
1088 cout << " qTotTPC "<< qTotTPC << endl;
1089 }
0cbdc86b 1090 //
1091 // glitch filter
1092 //
1093 digrow->GlitchFilter();
1094 //
1095 digrow->CompresBuffer(1,zerosup);
1096 digrow->CompresTrackBuffer(1);
1097 tree->Fill();
1098 if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
6fadf71b 1099 } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
1100
0cbdc86b 1101
1102 orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
1103 ogime = orl->GetLoader("TPCLoader");
1104 ogime->WriteDigits("OVERWRITE");
6fadf71b 1105
0cbdc86b 1106 //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
6fadf71b 1107
0cbdc86b 1108 delete digrow;
1109 for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
1110 delete []masks;
1111 delete []pdig;
1112 delete []ptr;
1113 delete []active;
1114 delete []digarr;
1115}