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