1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////////
18 // TRD MCM (Multi Chip Module) simulator //
20 ///////////////////////////////////////////////////////////////////////////////
26 New release on 2007/08/17
28 AliTRDmcmSim is now stably working and zero suppression function seems ok.
29 From now, the default version of raw data is set to 3 in AliTRDfeeParam.
31 The following internal parameters were abolished because it is useless and
37 GetCol member was modified accordingly.
39 New member function DumpData was prepared for diagnostics.
41 ZSMapping member function was debugged. It was causing crash due to
42 wrong indexing in 1 dimensional numbering. Also code was shaped up better.
46 /*Semi-final version of TRD raw data simulation code with zero suppression (ZS)
47 similar to TRD FEE. ZS is realized by the class group:
53 AliTRDfeeParam has been modified to have more parameters like raw data
54 production version and so on. AliTRDmcmSim is new class and this is the core
55 of MCM (PASA+TRAP) simulator. It has still very simple function and it will be
56 another project to improve this to make it closer to the reall FEE.
57 AliTRDrawData has been modified to use new class AliTRDmcmSim.
59 These modifications were tested on Aug. 02 HEAD version that code itself
60 compiles. I'm sure there must be still bugs and we need testing by as many as
61 possible persons now. Especially it seems HLT part is impacted by problems
62 because some parameters were moved from AliTRDrawData to AliTRDfeeParam (like
63 fRawVersion disappeared from AliTRDrawData).
65 In TRD definition, we have now 4 raw data versions.
67 0 very old offline version (by Bogdan)
68 1 test version (no zero suppression)
69 2 final version (no zero suppression)
70 3 test version (with zero suppression)
72 The default is still set to 2 in AliTRDfeeParam::fgkRAWversion and it uses
73 previously existing codes. If you set this to 3, AliTRDrawData changes behavior
74 to use AliTRDmcmSim with ZS.
76 Plan is after we make sure it works stably, we delete AliTRDmcm which is obsolete.
77 However it still take time because tracklet part is not yet touched.
78 The default raw version is 2.
83 // if no histo is drawn, these are obsolete
87 // only needed if I/O of tracklets is activated
99 #include "AliTRDmcmSim.h"
100 #include "AliTRDfeeParam.h"
101 #include "AliTRDSimParam.h"
102 #include "AliTRDgeometry.h"
103 #include "AliTRDcalibDB.h"
104 #include "AliTRDdigitsManager.h"
105 #include "AliTRDarrayADC.h"
106 // additional for new tail filter and/or tracklet
107 #include "AliTRDtrapAlu.h"
108 #include "AliTRDpadPlane.h"
109 #include "AliTRDtrackletMCM.h"
112 #include "AliLoader.h"
114 ClassImp(AliTRDmcmSim)
116 //_____________________________________________________________________________
117 AliTRDmcmSim::AliTRDmcmSim() :TObject()
118 ,fInitialized(kFALSE)
143 // AliTRDmcmSim default constructor
146 // By default, nothing is initialized.
147 // It is necessary to issue Init before use.
150 //_____________________________________________________________________________
151 AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m)
153 ,fInitialized(kFALSE)
179 // AliTRDmcmSim copy constructor
182 // By default, nothing is initialized.
183 // It is necessary to issue Init before use.
186 //_____________________________________________________________________________
187 AliTRDmcmSim::~AliTRDmcmSim()
190 // AliTRDmcmSim destructor
193 if( fADCR != NULL ) {
194 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
195 delete [] fADCR[iadc];
196 delete [] fADCF[iadc];
197 delete [] fADCT[iadc];
198 delete [] fZSM [iadc];
216 //_____________________________________________________________________________
217 AliTRDmcmSim &AliTRDmcmSim::operator=(const AliTRDmcmSim &m)
220 // Assignment operator
224 ((AliTRDmcmSim &) m).Copy(*this);
230 //_____________________________________________________________________________
231 void AliTRDmcmSim::Copy(TObject &m) const
236 ((AliTRDmcmSim &) m).fNextEvent = 0; //new
237 ((AliTRDmcmSim &) m).fMaxTracklets = 0; //new
238 ((AliTRDmcmSim &) m).fInitialized = 0;
239 ((AliTRDmcmSim &) m).fChaId = 0;
240 ((AliTRDmcmSim &) m).fSector = 0;
241 ((AliTRDmcmSim &) m).fStack = 0;
242 ((AliTRDmcmSim &) m).fLayer = 0;
243 ((AliTRDmcmSim &) m).fRobPos = 0;
244 ((AliTRDmcmSim &) m).fMcmPos = 0;
245 ((AliTRDmcmSim &) m).fNADC = 0;
246 ((AliTRDmcmSim &) m).fNTimeBin = 0;
247 ((AliTRDmcmSim &) m).fRow = 0;
248 ((AliTRDmcmSim &) m).fADCR = 0;
249 ((AliTRDmcmSim &) m).fADCF = 0;
250 ((AliTRDmcmSim &) m).fADCT = 0; //new
251 ((AliTRDmcmSim &) m).fPosLUT = 0; //new
252 ((AliTRDmcmSim &) m).fMCMT = 0; //new
253 ((AliTRDmcmSim &) m).fZSM = 0;
254 ((AliTRDmcmSim &) m).fZSM1Dim = 0;
255 ((AliTRDmcmSim &) m).fFeeParam = 0;
256 ((AliTRDmcmSim &) m).fSimParam = 0;
257 ((AliTRDmcmSim &) m).fCal = 0;
258 ((AliTRDmcmSim &) m).fGeo = 0;
262 //_____________________________________________________________________________
264 //void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos )
265 void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos, Bool_t newEvent = kFALSE ) // only for readout tree (new event)
268 // Initialize the class with new geometry information
269 // fADC array will be reused with filled by zero
273 fFeeParam = AliTRDfeeParam::Instance();
274 fSimParam = AliTRDSimParam::Instance();
275 fCal = AliTRDcalibDB::Instance();
276 fGeo = new AliTRDgeometry();
278 fSector = fGeo->GetSector( fChaId );
279 fStack = fGeo->GetStack( fChaId );
280 fLayer = fGeo->GetLayer( fChaId );
283 fNADC = fFeeParam->GetNadcMcm();
284 fNTimeBin = fCal->GetNumberOfTimeBins();
285 fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
287 fMaxTracklets = fFeeParam->GetMaxNrOfTracklets();
292 if (newEvent == kTRUE) {
298 // Allocate ADC data memory if not yet done
299 if( fADCR == NULL ) {
300 fADCR = new Int_t *[fNADC];
301 fADCF = new Int_t *[fNADC];
302 fADCT = new Int_t *[fNADC]; //new
303 fZSM = new Int_t *[fNADC];
304 fZSM1Dim = new Int_t [fNADC];
305 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
306 fADCR[iadc] = new Int_t[fNTimeBin];
307 fADCF[iadc] = new Int_t[fNTimeBin];
308 fADCT[iadc] = new Int_t[fNTimeBin]; //new
309 fZSM [iadc] = new Int_t[fNTimeBin];
313 // Initialize ADC data
314 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
315 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
318 fADCT[iadc][it] = -1; //new
319 fZSM [iadc][it] = 1; // Default unread = 1
321 fZSM1Dim[iadc] = 1; // Default unread = 1
325 fPosLUT = new Int_t[128];
326 for(Int_t i = 0; i<128; i++){
330 fMCMT = new UInt_t[fMaxTracklets];
331 for(Int_t i = 0; i < fMaxTracklets; i++) {
336 fInitialized = kTRUE;
339 //_____________________________________________________________________________
340 Bool_t AliTRDmcmSim::CheckInitialized()
343 // Check whether object is initialized
346 if( ! fInitialized ) {
347 AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
352 //_____________________________________________________________________________
355 void AliTRDmcmSim::SetPosLUT() {
356 Double_t iHi = (Double_t)fCal->GetPRFhi();
357 Double_t iLo = (Double_t)fCal->GetPRFlo();
358 Int_t nBin = fCal->GetPRFbin();
359 Int_t iOff = fLayer * nBin;
360 Int_t kNlayer = fGeo->Nlayer();
362 Float_t *sPRFsmp = new Float_t[nBin*kNlayer];
363 Double_t *sPRFlayer = new Double_t[nBin];
366 for(Int_t i = 0; i<nBin*kNlayer; i++){
368 //printf("%f\n",fCal->GetSampledPRF()[i]);
369 sPRFsmp[i] = fCal->GetSampledPRF()[i];
373 Double_t sWidth = (iHi-iLo)/((Double_t) nBin);
374 Int_t sPad = (Int_t) (1.0/sWidth);
376 // get the PRF for actual layer (interpolated to ibin data-points; 61 measured)
377 for(Int_t iBin = 0; iBin < nBin; iBin++){
378 sPRFlayer[iBin] = (Double_t)sPRFsmp[iOff+iBin];
381 Int_t bin0 = (Int_t)(-iLo / sWidth - 0.5); // bin-nr. for pad-position 0
383 Int_t bin1 = (Int_t)((Double_t)(0.5 - iLo) / sWidth - 0.5); // bin-nr. for pad-position 0.5
385 bin0 = bin0 + 1; //avoid negative values in aYest (start right of symmetry center)
386 while (bin0-sPad<0) {
389 while (bin1+sPad>=nBin) {
393 Double_t* aYest = new Double_t[bin1-bin0+1];
395 /*TH1F* hist1 = new TH1F("h1","yest(y)",128,0,0.5);
396 TH1F* hist2 = new TH1F("h2","y(yest)",128,0,0.5);
397 TH1F* hist3 = new TH1F("h3","y(yest)-yest",128,0,0.5);
398 TH1F* hist4 = new TH1F("h4","y(yest)-yest,discrete",128,0,0.5);
400 TCanvas *c1 = new TCanvas("c1","c1",800,1000);
402 TCanvas *c2 = new TCanvas("c2","c2",800,1000);
404 TCanvas *c3 = new TCanvas("c3","c3",800,1000);
406 TCanvas *c4 = new TCanvas("c4","c4",800,1000);
409 for(Int_t iBin = bin0; iBin <= bin1; iBin++){
410 aYest[iBin-bin0] = 0.5*(sPRFlayer[iBin-sPad] - sPRFlayer[iBin+sPad])/(sPRFlayer[iBin]); // estimated position from PRF; between 0 and 1
411 //Double_t position = ((Double_t)(iBin)+0.5)*sWidth+iLo;
412 // hist1->Fill(position,aYest[iBin-bin0]);
417 Double_t aY[128]; // reversed function
422 for(Int_t j = 0; j<128; j++) { // loop over all Yest; LUT has 128 entries;
423 Double_t yest = ((Double_t)j)/256;
426 while (yest>aYest[iBin] && iBin<(bin1-bin0)) {
429 if((iBin == bin1 - bin0)&&(yest>aYest[iBin])) {
430 aY[j] = 0.5; // yest too big
431 //hist2->Fill(yest,aY[j]);
435 Int_t bin_d = iBin + bin0 - 1;
436 Int_t bin_u = iBin + bin0;
437 Double_t y_d = ((Double_t)bin_d + 0.5)*sWidth + iLo; // lower y
438 Double_t y_u = ((Double_t)bin_u + 0.5)*sWidth + iLo; // upper y
439 Double_t yest_d = aYest[iBin-1]; // lower estimated y
440 Double_t yest_u = aYest[iBin]; // upper estimated y
442 aY[j] = ((yest-yest_d)/(yest_u-yest_d))*(y_u-y_d) + y_d;
443 //hist2->Fill(yest,aY[j]);
446 aY[j] = aY[j] - yest;
447 //hist3->Fill(yest,aY[j]);
449 a.AssignDouble(aY[j]);
451 fPosLUT[j] = a.GetValue(); // 1+8Bit value;128 entries;LUT is steered by abs(Q(i+1)-Q(i-1))/Q(i)=COG and gives the correction to COG/2
452 //hist4->Fill(yest,fPosLUT[j]);
465 //_____________________________________________________________________________
466 Int_t* AliTRDmcmSim::GetPosLUT(){
472 void AliTRDmcmSim::SetData( Int_t iadc, Int_t *adc )
475 // Store ADC data into array of raw data
478 if( !CheckInitialized() ) return;
480 if( iadc < 0 || iadc >= fNADC ) {
481 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
485 for( int it = 0 ; it < fNTimeBin ; it++ ) {
486 fADCR[iadc][it] = (Int_t)(adc[it]);
490 //_____________________________________________________________________________
491 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
494 // Store ADC data into array of raw data
497 if( !CheckInitialized() ) return;
499 if( iadc < 0 || iadc >= fNADC ) {
500 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
504 fADCR[iadc][it] = adc;
507 //_____________________________________________________________________________
508 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
511 // Store ADC data into array of raw data
514 if( !CheckInitialized() ) return;
516 if( iadc < 0 || iadc >= fNADC ) {
517 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
521 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
522 fADCR[iadc][it] = fSimParam->GetADCbaseline();
526 //_____________________________________________________________________________
527 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
530 // Return column id of the pad for the given ADC channel
533 if( !CheckInitialized() ) return -1;
535 return fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
538 //_____________________________________________________________________________
539 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize )
542 // Produce raw data stream from this MCM and put in buf
543 // Returns number of words filled, or negative value
544 // with -1 * number of overflowed words
549 Int_t nw = 0; // Number of written words
550 Int_t of = 0; // Number of overflowed words
551 Int_t rawVer = fFeeParam->GetRAWversion();
554 if( !CheckInitialized() ) return 0;
556 if( fFeeParam->GetRAWstoreRaw() ) {
562 // Produce MCM header
563 x = ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
574 for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
575 if( fZSM1Dim[iAdc] == 0 ) { // 0 means not suppressed
587 // Produce ADC data. 3 timebins are packed into one 32 bits word
588 // In this version, different ADC channel will NOT share the same word
590 UInt_t aa=0, a1=0, a2=0, a3=0;
592 for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
593 if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // suppressed
594 aa = !(iAdc & 1) + 2;
595 for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
596 a1 = ((iT ) < fNTimeBin ) ? adc[iAdc][iT ] : 0;
597 a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
598 a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
599 x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
609 if( of != 0 ) return -of; else return nw;
612 //_____________________________________________________________________________
613 Int_t AliTRDmcmSim::ProduceRawStreamV2( UInt_t *buf, Int_t maxSize )
616 // Produce raw data stream from this MCM and put in buf
617 // Returns number of words filled, or negative value
618 // with -1 * number of overflowed words
623 Int_t nw = 0; // Number of written words
624 Int_t of = 0; // Number of overflowed words
625 Int_t rawVer = fFeeParam->GetRAWversion();
627 Int_t nActiveADC = 0; // number of activated ADC bits in a word
629 if( !CheckInitialized() ) return 0;
631 if( fFeeParam->GetRAWstoreRaw() ) {
637 // Produce MCM header
638 x = (1<<31) | ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
641 //printf("\nMCM header: %X ",x);
647 // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
648 // n : unused , c : ADC count, m : selected ADCs
651 for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
652 if( fZSM1Dim[iAdc] == 0 ) { // 0 means not suppressed
653 x = x | (1 << (iAdc+4) ); // last 4 digit reserved for 1100=0xc
654 nActiveADC++; // number of 1 in mmm....m
657 x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC; // nn = 01, ccccc are inverted, 0xc=1100
658 //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x );
662 //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
669 // Produce ADC data. 3 timebins are packed into one 32 bits word
670 // In this version, different ADC channel will NOT share the same word
672 UInt_t aa=0, a1=0, a2=0, a3=0;
674 for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
675 if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
676 aa = !(iAdc & 1) + 2;
677 for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
678 a1 = ((iT ) < fNTimeBin ) ? adc[iAdc][iT ] : 0;
679 a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
680 a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
681 x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
692 if( of != 0 ) return -of; else return nw;
695 //_____________________________________________________________________________
696 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
699 // Produce tracklet data stream from this MCM and put in buf
700 // Returns number of words filled, or negative value
701 // with -1 * number of overflowed words
705 Int_t nw = 0; // Number of written words
706 Int_t of = 0; // Number of overflowed words
708 if( !CheckInitialized() ) return 0;
710 // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM
711 // fMCMT is filled continuously until no more tracklet words available
714 while ( (wd < fMaxTracklets) && (fMCMT[wd] > 0) ){
725 if( of != 0 ) return -of; else return nw;
729 //_____________________________________________________________________________
730 void AliTRDmcmSim::Filter()
733 // Apply digital filter
736 if( !CheckInitialized() ) return;
738 // Initialize filtered data array with raw data
739 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
740 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
741 fADCF[iadc][it] = fADCR[iadc][it];
745 // Then apply fileters one by one to filtered data array
746 if( fFeeParam->IsPFon() ) FilterPedestal();
747 if( fFeeParam->IsGFon() ) FilterGain();
748 if( fFeeParam->IsTFon() ) FilterTail();
751 //_____________________________________________________________________________
752 void AliTRDmcmSim::FilterPedestal()
757 // Apply pedestal filter
760 Int_t ap = fSimParam->GetADCbaseline(); // ADC instrinsic pedestal
761 Int_t ep = fFeeParam->GetPFeffectPedestal(); // effective pedestal
762 //Int_t tc = fFeeParam->GetPFtimeConstant(); // this makes no sense yet
764 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
765 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
766 fADCF[iadc][it] = fADCF[iadc][it] - ap + ep;
771 //_____________________________________________________________________________
772 void AliTRDmcmSim::FilterGain()
775 // Apply gain filter (not implemented)
776 // Later it will be implemented because gain digital filiter will
777 // increase noise level.
782 //_____________________________________________________________________________
783 void AliTRDmcmSim::FilterTail()
786 // Apply exponential tail filter (Bogdan's version)
789 Double_t *dtarg = new Double_t[fNTimeBin];
790 Int_t *itarg = new Int_t[fNTimeBin];
791 Int_t nexp = fFeeParam->GetTFnExp();
792 Int_t tftype = fFeeParam->GetTFtype();
796 case 0: // Exponential Filter Analog Bogdan
797 for (Int_t iCol = 0; iCol < fNADC; iCol++) {
798 FilterSimDeConvExpA( fADCF[iCol], dtarg, fNTimeBin, nexp);
799 for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
800 fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
805 case 1: // Exponential filter digital Bogdan
806 for (Int_t iCol = 0; iCol < fNADC; iCol++) {
807 FilterSimDeConvExpD( fADCF[iCol], itarg, fNTimeBin, nexp);
808 for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
809 fADCF[iCol][iTime] = itarg[iTime];
814 case 2: // Exponential filter Marian special
815 for (Int_t iCol = 0; iCol < fNADC; iCol++) {
816 FilterSimDeConvExpMI( fADCF[iCol], dtarg, fNTimeBin);
817 for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
818 fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
824 case 3: // Exponential filter using AliTRDtrapAlu class
825 for (Int_t iCol = 0; iCol < fNADC; iCol++) {
826 FilterSimDeConvExpEl( fADCF[iCol], itarg, fNTimeBin, nexp);
827 for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
828 fADCF[iCol][iTime] = itarg[iTime]>>2; // to be used for raw-data
829 fADCT[iCol][iTime] = itarg[iTime]; // 12bits; to be used for tracklet; tracklet will have own container;
836 AliError(Form("Invalid filter type %d ! \n", tftype ));
844 //_____________________________________________________________________________
845 void AliTRDmcmSim::ZSMapping()
848 // Zero Suppression Mapping implemented in TRAP chip
850 // See detail TRAP manual "Data Indication" section:
851 // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
854 Int_t eBIS = fFeeParam->GetEBsglIndThr(); // TRAP default = 0x4 (Tis=4)
855 Int_t eBIT = fFeeParam->GetEBsumIndThr(); // TRAP default = 0x28 (Tit=40)
856 Int_t eBIL = fFeeParam->GetEBindLUT(); // TRAP default = 0xf0
857 // (lookup table accept (I2,I1,I0)=(111)
858 // or (110) or (101) or (100))
859 Int_t eBIN = fFeeParam->GetEBignoreNeighbour(); // TRAP default = 1 (no neighbor sensitivity)
860 Int_t ep = AliTRDfeeParam::GetPFeffectPedestal();
862 if( !CheckInitialized() ) return;
864 for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
865 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
867 // Get ADC data currently in filter buffer
868 Int_t ap = fADCF[iadc-1][it] - ep; // previous
869 Int_t ac = fADCF[iadc ][it] - ep; // current
870 Int_t an = fADCF[iadc+1][it] - ep; // next
872 // evaluate three conditions
873 Int_t i0 = ( ac >= ap && ac >= an ) ? 0 : 1; // peak center detection
874 Int_t i1 = ( ap + ac + an > eBIT ) ? 0 : 1; // cluster
875 Int_t i2 = ( ac > eBIS ) ? 0 : 1; // absolute large peak
877 Int_t i = i2 * 4 + i1 * 2 + i0; // Bit position in lookup table
878 Int_t d = (eBIL >> i) & 1; // Looking up (here d=0 means true
879 // and d=1 means false according to TRAP manual)
882 if( eBIN == 0 ) { // turn on neighboring ADCs
883 fZSM[iadc-1][it] &= d;
884 fZSM[iadc+1][it] &= d;
890 // do 1 dim projection
891 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
892 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
893 fZSM1Dim[iadc] &= fZSM[iadc][it];
899 //_____________________________________________________________________________
900 void AliTRDmcmSim::DumpData( char *f, char *target )
903 // Dump data stored (for debugging).
904 // target should contain one or multiple of the following characters
906 // F for filtered data
907 // Z for zero suppression map
909 // other characters are simply ignored
912 UInt_t tempbuf[1024];
914 if( !CheckInitialized() ) return;
916 std::ofstream of( f, std::ios::out | std::ios::app );
917 of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
918 fChaId, fSector, fStack, fLayer, fRobPos, fMcmPos );
920 for( int t=0 ; target[t] != 0 ; t++ ) {
921 switch( target[t] ) {
924 of << Form("fADCR (raw ADC data)\n");
925 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
926 of << Form(" ADC %02d: ", iadc);
927 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
928 of << Form("% 4d", fADCR[iadc][it]);
935 of << Form("fADCF (filtered ADC data)\n");
936 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
937 of << Form(" ADC %02d: ", iadc);
938 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
939 of << Form("% 4d", fADCF[iadc][it]);
946 of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
947 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
948 of << Form(" ADC %02d: ", iadc);
949 if( fZSM1Dim[iadc] == 0 ) { of << " R " ; } else { of << " . "; } // R:read .:suppressed
950 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
951 if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
958 Int_t s = ProduceRawStream( tempbuf, 1024 );
959 of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
960 of << Form(" address data\n");
961 for( int i = 0 ; i < s ; i++ ) {
962 of << Form(" %04x %08x\n", i, tempbuf[i]);
968 //_____________________________________________________________________________
969 void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target
970 , Int_t n, Int_t nexp)
973 // Exponential filter "analog"
974 // source will not be changed
979 Double_t reminder[2];
983 Double_t coefficients[2];
985 // Initialize (coefficient = alpha, rates = lambda)
986 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
988 Double_t r1 = (Double_t)fFeeParam->GetTFr1();
989 Double_t r2 = (Double_t)fFeeParam->GetTFr2();
990 Double_t c1 = (Double_t)fFeeParam->GetTFc1();
991 Double_t c2 = (Double_t)fFeeParam->GetTFc2();
993 coefficients[0] = c1;
994 coefficients[1] = c2;
997 rates[0] = TMath::Exp(-dt/(r1));
998 rates[1] = TMath::Exp(-dt/(r2));
1000 // Attention: computation order is important
1002 for (k = 0; k < nexp; k++) {
1006 for (i = 0; i < n; i++) {
1008 result = ((Double_t)source[i] - correction); // no rescaling
1011 for (k = 0; k < nexp; k++) {
1012 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1016 for (k = 0; k < nexp; k++) {
1017 correction += reminder[k];
1022 //_____________________________________________________________________________
1023 void AliTRDmcmSim::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n
1027 // Exponential filter "digital"
1028 // source will not be changed
1038 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
1039 // initialize (coefficient = alpha, rates = lambda)
1042 Double_t r1 = (Double_t)fFeeParam->GetTFr1();
1043 Double_t r2 = (Double_t)fFeeParam->GetTFr2();
1044 Double_t c1 = (Double_t)fFeeParam->GetTFc1();
1045 Double_t c2 = (Double_t)fFeeParam->GetTFc2();
1047 Int_t fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
1048 Int_t fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
1049 Int_t iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; // 9 bit paramter + fixed bits
1050 Int_t iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; // 9 bit paramter + fixed bits
1053 fAlphaL = (Int_t) (c1 * 2048.0);
1054 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
1057 fAlphaL = (Int_t) (c1 * 2048.0);
1058 fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
1059 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
1060 iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400; // 10 bit paramter + fixed bits
1063 Double_t iAl = iAlphaL / 2048.0; // alpha L: correspondence to floating point numbers
1064 Double_t iAs = iAlphaS / 2048.0; // alpha S: correspondence to floating point numbers
1065 Double_t iLl = iLambdaL / 2048.0; // lambda L: correspondence to floating point numbers
1066 Double_t iLs = iLambdaS / 2048.0; // lambda S: correspondence to floating point numbers
1074 Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2;
1076 Double_t xi = 1 - (iLl*iAs + iLs*iAl); // Calculation of equilibrium values of the
1077 rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
1078 rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
1080 // further initialization
1081 if ((rem1 + rem2) > 0x0FFF) {
1082 correction = 0x0FFF;
1085 correction = (rem1 + rem2) & 0x0FFF;
1088 fTailPed = iFactor - correction;
1090 for (i = 0; i < n; i++) {
1092 result = (source[i] - correction);
1093 if (result < 0) { // Too much undershoot
1099 h1 = (rem1 + ((iAlphaL * result) >> 11));
1107 h2 = (rem2 + ((iAlphaS * result) >> 11));
1115 rem1 = (iLambdaL * h1 ) >> 11;
1116 rem2 = (iLambdaS * h2 ) >> 11;
1118 if ((rem1 + rem2) > 0x0FFF) {
1119 correction = 0x0FFF;
1122 correction = (rem1 + rem2) & 0x0FFF;
1129 //_____________________________________________________________________________
1130 void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target
1134 // Exponential filter (M. Ivanov)
1135 // source will not be changed
1143 for (i = 0; i < n; i++) {
1144 sig1[i] = (Double_t)source[i];
1148 Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt;
1149 Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt;
1151 FilterSimTailMakerSpline( sig1, sig2, lambda0, n);
1152 FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n);
1154 for (i = 0; i < n; i++) {
1155 target[i] = sig3[i];
1160 //______________________________________________________________________________
1161 void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout
1162 , Double_t lambda, Int_t n)
1165 // Special filter (M. Ivanov)
1169 Double_t l = TMath::Exp(-lambda*0.5);
1173 // Initialize in[] and out[] goes 0 ... 2*n+19
1174 for (i = 0; i < n*2+20; i++) {
1180 in[1] = (ampin[0] + ampin[1]) * 0.5;
1182 // Add charge to the end
1183 for (i = 0; i < 22; i++) {
1184 in[2*(n-1)+i] = ampin[n-1]; // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
1187 // Use arithmetic mean
1188 for (i = 1; i < n-1; i++) {
1189 in[2*i] = ampin[i]; // in[] goes 2, 3, ... , 2*n-4, 2*n-3
1190 in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
1196 for (i = 2*n; i >= 0; i--) {
1197 out[i] = in[i] + temp;
1198 temp = l*(temp+in[i]);
1201 for (i = 0; i < n; i++){
1202 //ampout[i] = out[2*i+1]; // org
1203 ampout[i] = out[2*i];
1208 //______________________________________________________________________________
1209 void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout
1210 , Double_t norm, Double_t lambda
1214 // Special filter (M. Ivanov)
1219 Double_t l = TMath::Exp(-lambda*0.5);
1220 Double_t k = l*(1.0 - norm*lambda*0.5);
1224 // Initialize in[] and out[] goes 0 ... 2*n+19
1225 for (i = 0; i < n*2+20; i++) {
1231 in[1] = (ampin[0]+ampin[1])*0.5;
1233 // Add charge to the end
1234 for (i =-2; i < 22; i++) {
1235 // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
1236 in[2*(n-1)+i] = ampin[n-1];
1239 for (i = 1; i < n-2; i++) {
1240 // in[] goes 2, 3, ... , 2*n-6, 2*n-5
1242 in[2*i+1] = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
1243 //in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.0;
1249 for (i = 1; i <= 2*n; i++) {
1250 out[i] = in[i] + (k-l)*temp;
1251 temp = in[i] + k *temp;
1254 for (i = 0; i < n; i++) {
1255 //ampout[i] = out[2*i+1]; // org
1256 //ampout[i] = TMath::Max(out[2*i+1],0.0); // org
1257 ampout[i] = TMath::Max(out[2*i],0.0);
1262 //_____________________________________________________________________________________
1263 //the following filter uses AliTRDtrapAlu-class
1265 void AliTRDmcmSim::FilterSimDeConvExpEl(Int_t *source, Int_t *target, Int_t n, Int_t nexp) {
1266 //static Int_t count = 0;
1269 Double_t r1 = (Double_t)fFeeParam->GetTFr1();
1270 Double_t r2 = (Double_t)fFeeParam->GetTFr2();
1271 Double_t c1 = (Double_t)fFeeParam->GetTFc1();
1272 Double_t c2 = (Double_t)fFeeParam->GetTFc2();
1276 //it is assumed that r1,r2,c1,c2 are given such, that the configuration values are in the ranges according to TRAP-manual
1277 //parameters need to be adjusted
1278 AliTRDtrapAlu lambdaL;
1279 AliTRDtrapAlu lambdaS;
1280 AliTRDtrapAlu alphaL;
1281 AliTRDtrapAlu alphaS;
1283 AliTRDtrapAlu correction;
1284 AliTRDtrapAlu result;
1288 AliTRDtrapAlu bSource;
1297 lambdaL.AssignDouble(TMath::Exp(-dt/r1));
1298 lambdaS.AssignDouble(TMath::Exp(-dt/r2));
1299 alphaL.AssignDouble(c1); // in AliTRDfeeParam the number of exponentials is set and also the according time constants
1300 alphaS.AssignDouble(c2); // later it should be: alphaS=1-alphaL
1302 //data is enlarged to 12 bits, including 2 bits after the comma; class AliTRDtrapAlu is used to handle arithmetics correctly
1303 correction.Init(10,2);
1309 for(Int_t i = 0; i < n; i++) {
1310 bSource.AssignInt(source[i]);
1311 result = bSource - correction; // subtraction can produce an underflow
1312 if(result.GetSign() == kTRUE) {
1313 result.AssignInt(0);
1316 //target[i] = result.GetValuePre(); // later, target and source should become AliTRDtrapAlu,too in order to simulate the 10+2Bits through the filter properly
1318 target[i] = result.GetValue(); // 12 bit-value; to get the corresponding integer value, target must be shifted: target>>2
1320 //printf("target-Wert zur Zeit %d : %d",i,target[i]);
1323 bufL = bufL + (result * alphaL);
1324 bufL = bufL * lambdaL;
1326 bufS = bufS + (result * alphaS);
1327 bufS = bufS * lambdaS; // eventually this should look like:
1328 // bufS = (bufS + (result - result * alphaL)) * lambdaS // alphaS=1-alphaL; then alphaS-variable is not needed any more
1330 correction = bufL + bufS; //check for overflow intrinsic; if overflowed, correction is set to 0x03FF
1342 //__________________________________________________________________________________
1345 // in order to use the Tracklets, please first
1346 // -- set AliTRDfeeParam::fgkTracklet to kTRUE, in order to switch on Tracklet-calculation
1347 // -- set AliTRDfeeParam::fgkTFtype to 3, in order to use the new tail cancellation filter
1348 // currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3
1349 // -- set AliTRDfeeParam::fgkMCTrackletOutput to kTRUE, if you want to use the Tracklet output container with information about the Tracklet position (MCM, channel number)
1351 // The code is designed such that the less possible calculations with AliTRDtrapAlu class-objects are performed; whenever possible calculations are done with doubles or integers and the results are transformed into the right format
1353 void AliTRDmcmSim::Tracklet(){
1354 // tracklet calculation
1355 // if you use this code after a simulation, please make sure the same filter-settings as in the simulation are set in AliTRDfeeParam
1357 if(!CheckInitialized()){ return; }
1359 Bool_t filtered = kTRUE;
1365 if(fADCT[0][0]==-1){ // check if filter was applied
1367 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1368 for( Int_t iT = 0 ; iT < fNTimeBin ; iT++ ) {
1369 data.AssignInt(fADCR[iadc][iT]);
1370 fADCT[iadc][iT] = data.GetValue(); // all incoming values are positive 10+2 bit values; if el.filter was called, this is done correctly
1376 // the online ordering of mcm's is reverse to the TRAP-manual-ordering! reverse fADCT (to be consistent to TRAP), then do all calculations
1378 Int_t** rev0 = new Int_t *[fNADC];
1379 Int_t** rev1 = new Int_t *[fNADC];
1381 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1382 rev0[iadc] = new Int_t[fNTimeBin];
1383 rev1[iadc] = new Int_t[fNTimeBin];
1384 for( Int_t iT = 0; iT < fNTimeBin; iT++) {
1385 if( iadc <= fNADC-iadc-1 ) {
1386 rev0[iadc][iT] = fADCT[fNADC-iadc-1][iT];
1387 rev1[iadc][iT] = fADCT[iadc][iT];
1388 fADCT[iadc][iT] = rev0[iadc][iT];
1391 rev0[iadc][iT] = rev1[fNADC-iadc-1][iT];
1392 fADCT[iadc][iT] = rev0[iadc][iT];
1396 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1397 delete[] rev0[iadc];
1398 delete[] rev1[iadc];
1407 // get the filtered pedestal; supports only electronic tail-cancellation filter
1408 AliTRDtrapAlu filPed;
1410 Int_t *ieffped = new Int_t[fNTimeBin];
1411 for(Int_t iT = 0; iT < fNTimeBin; iT++){
1415 if( filtered == kTRUE ) {
1416 if( fFeeParam->IsPFon() ){
1417 ep = fFeeParam->GetPFeffectPedestal();
1419 Int_t nexp = fFeeParam->GetTFnExp();
1420 Int_t *isource = new Int_t[fNTimeBin];
1422 filPed.AssignInt(ep);
1423 Int_t epf = filPed.GetValue();
1424 for(Int_t iT = 0; iT < fNTimeBin; iT++){
1429 if( fFeeParam->IsTFon() ) {
1430 FilterSimDeConvExpEl( isource, ieffped, fNTimeBin, nexp);
1436 //the following values should go to AliTRDfeeParam once they are defined; then they have to be read in properly
1437 //naming follows conventions in TRAP-manual
1440 Bool_t bVBY = kTRUE; // cluster-verification bypass
1442 Double_t cQTParam = 0; // cluster quality threshold; granularity 2^-10; range: 0<=cQT/2^-10<=2^-4 - 2^-10
1443 AliTRDtrapAlu cQTAlu;
1444 cQTAlu.Init(1,10,0,63);
1445 cQTAlu.AssignDouble(cQTParam);
1446 Int_t cQT = cQTAlu.GetValue();
1449 Int_t tFS = fFeeParam->GetLinearFitStart(); // linear fit start
1450 Int_t tFE = fFeeParam->GetLinearFitEnd(); // linear fit stop
1452 // charge accumulators
1453 Int_t tQS0 = fFeeParam->GetQacc0Start(); // start-time for charge-accumulator 0
1454 Int_t tQE0 = fFeeParam->GetQacc0End(); // stop-time for charge-accumulator 0
1455 Int_t tQS1 = fFeeParam->GetQacc1Start(); // start-time for charge-accumulator 1
1456 Int_t tQE1 = fFeeParam->GetQacc1End(); // stop-time for charge-accumulator 1
1457 // values set such that tQS0=tFS; tQE0=tQS1-1; tFE=tQE1; want to do (QS0+QS1)/N
1459 Double_t cTHParam = (Double_t)fFeeParam->GetMinClusterCharge(); // cluster charge threshold
1460 AliTRDtrapAlu cTHAlu;
1462 cTHAlu.AssignDouble(cTHParam);
1463 Int_t cTH = cTHAlu.GetValue(); // cTH used for comparison
1471 List_t selection[7]; // list with 7 elements
1472 List_t *list = NULL;
1473 List_t *listLeft = NULL;
1475 Int_t* qsum = new Int_t[fNADC];
1478 AliTRDtrapAlu qsumAlu;
1479 qsumAlu.Init(12,2); // charge sum will be 12+2 bits
1480 AliTRDtrapAlu dCOGAlu;
1481 dCOGAlu.Init(1,7,0,127); // COG will be 1+7 Bits; maximum 1 - 2^-7 for LUT
1482 AliTRDtrapAlu yrawAlu;
1483 yrawAlu.Init(1,8,-1,255);
1485 yAlu.Init(1,16,-1,0xFF00); // only first 8 past-comma bits filled;additional 8 bits for accuracy;maximum 1 - 2^-8; sign is given by + or -
1487 xAlu.Init(5,8); // 8 past-comma bits because value will be added/multiplied to another value with this accuracy
1488 AliTRDtrapAlu xxAlu;
1490 AliTRDtrapAlu yyAlu;
1491 yyAlu.Init(1,16,0,0xFFFF); // maximum is 2^16-1; 16Bit for past-commas
1492 AliTRDtrapAlu xyAlu;
1496 AliTRDtrapAlu XXAlu;
1499 YAlu.Init(5,8); // 14 bit, 1 is sign-bit; therefore only 13 bit
1500 AliTRDtrapAlu YYAlu;
1502 AliTRDtrapAlu XYAlu;
1503 XYAlu.Init(8,8); // 17 bit, 1 is sign-bit; therefore only 16 bit
1504 AliTRDtrapAlu qtruncAlu;
1505 qtruncAlu.Init(12,0);
1506 AliTRDtrapAlu QT0Alu;
1508 AliTRDtrapAlu QT1Alu;
1511 AliTRDtrapAlu oneAlu;
1515 AliTRDtrapAlu inverseNAlu;
1516 inverseNAlu.Init(1,8); // simulates the LUT for 1/N
1517 AliTRDtrapAlu MeanChargeAlu; // mean charge in ADC counts
1518 MeanChargeAlu.Init(8,0);
1519 AliTRDtrapAlu TotalChargeAlu;
1520 TotalChargeAlu.Init(17,8);
1521 //nr of post comma bits should be the same for inverseN and TotalCharge
1524 SetPosLUT(); // initialize the position correction LUT for this MCM;
1527 // fit-sums; remapping!; 0,1,2->0; 1,2,3->1; ... 18,19,20->18
1528 Int_t *X = new Int_t[fNADC-2];
1529 Int_t *XX = new Int_t[fNADC-2];
1530 Int_t *Y = new Int_t[fNADC-2];
1531 Int_t *YY = new Int_t[fNADC-2];
1532 Int_t *XY = new Int_t[fNADC-2];
1533 Int_t *N = new Int_t[fNADC-2];
1534 Int_t *QT0 = new Int_t[fNADC-2]; // accumulated charge
1535 Int_t *QT1 = new Int_t[fNADC-2]; // accumulated charge
1537 for (Int_t iCol = 0; iCol < fNADC-2; iCol++) {
1539 // initialize fit-sums
1551 filPed.Init(7,2); // convert filtered pedestal into 7+2Bits
1553 for(Int_t iT = 0; iT < fNTimeBin; iT++){
1555 if(iT<tFS || iT>=tFE) continue; // linear fit yes/no?
1558 Int_t portChannel[4] = {-1,-1,-1,-1};
1559 Int_t clusterCharge[4] = {0,0,0,0};
1560 Int_t leftCharge[4] = {0,0,0,0};
1561 Int_t centerCharge[4] = {0,0,0,0};
1562 Int_t rightCharge[4] = {0,0,0,0};
1566 filPed.AssignFormatted(ieffped[iT]); // no size-checking when using AssignFormatted; ieffped>=0
1567 filPed = filPed; // this checks the size
1569 ieffped[iT] = filPed.GetValue();
1571 for(Int_t i = 0; i<7; i++){
1572 selection[i].next = NULL;
1573 selection[i].iadc = -1; // value of -1: invalid adc
1574 selection[i].value = 0;
1577 // selection[0] is starting list-element; just for pointing
1579 // loop over inner adc's
1580 for (Int_t iCol = 1; iCol < fNADC-1; iCol++) {
1582 Int_t left = fADCT[iCol-1][iT];
1583 Int_t center = fADCT[iCol][iT];
1584 Int_t right = fADCT[iCol+1][iT];
1586 Int_t sum = left + center + right; // cluster charge sum
1587 qsumAlu.AssignFormatted(sum);
1588 qsumAlu = qsumAlu; // size-checking; redundant
1590 qsum[iCol] = qsumAlu.GetValue();
1592 //hit detection and masking
1595 if(qsum[iCol]>=(cTH + 3*ieffped[iT])){ // effective pedestal of all three channels must be added to cTH(+20); this is not parallel to TRAP manual; maybe cTH has to be adjusted in fFeeParam; therefore channels are not yet reduced by their pedestal
1596 mark |= 1; // marker
1605 // get selection of 6 adc's and sort,starting with greatest values
1607 //read three from right side and sort (primitive sorting algorithm)
1608 Int_t i = 0; // adc number
1609 Int_t j = 1; // selection number
1610 while(i<fNADC-2 && j<=3){
1612 if( ((mark>>(i-1)) & 1) == 1) {
1613 selection[j].iadc = fNADC-1-i;
1614 selection[j].value = qsum[fNADC-1-i]>>6; // for hit-selection only the first 8 out of the 14 Bits are used for comparison
1616 // insert into sorted list
1617 listLeft = &selection[0];
1618 list = listLeft->next;
1621 while((list->next != NULL) && (selection[j].value <= list->value)){
1626 if(selection[j].value<=list->value){
1627 selection[j].next = list->next;
1628 list->next = &selection[j];
1631 listLeft->next = &selection[j];
1632 selection[j].next = list;
1636 listLeft->next = &selection[j];
1637 selection[j].next = list;
1645 // read three from left side
1647 while(k>i && j<=6) {
1648 if( ((mark>>(k-1)) & 1) == 1) {
1649 selection[j].iadc = fNADC-1-k;
1650 selection[j].value = qsum[fNADC-1-k]>>6;
1652 listLeft = &selection[0];
1653 list = listLeft->next;
1656 while((list->next != NULL) && (selection[j].value <= list->value)){
1661 if(selection[j].value<=list->value){
1662 selection[j].next = list->next;
1663 list->next = &selection[j];
1666 listLeft->next = &selection[j];
1667 selection[j].next = list;
1671 listLeft->next = &selection[j];
1672 selection[j].next = list;
1680 // get the four with greatest charge-sum
1681 list = &selection[0];
1682 for(i = 0; i<4; i++){
1683 if(list->next == NULL) continue;
1685 if(list->iadc == -1) continue;
1686 Int_t adc = list->iadc; // channel number with selected hit
1688 // the following arrays contain the four chosen channels in 1 time-bin
1689 portChannel[i] = adc;
1690 clusterCharge[i] = qsum[adc];
1691 leftCharge[i] = fADCT[adc-1][iT] - ieffped[iT]; // reduce by filtered pedestal (pedestal is part of the signal)
1692 centerCharge[i] = fADCT[adc][iT] - ieffped[iT];
1693 rightCharge[i] = fADCT[adc+1][iT] - ieffped[iT];
1698 // cluster verification
1700 for(i = 0; i<4; i++){
1701 Int_t lr = leftCharge[i]*rightCharge[i]*1024;
1702 Int_t cc = centerCharge[i]*centerCharge[i]*cQT;
1704 portChannel[i] = -1; // set to invalid address
1705 clusterCharge[i] = 0;
1710 // fit-sums of valid channels
1711 // local hit position
1712 for(i = 0; i<4; i++){
1713 if (centerCharge[i] == 0) {
1714 portChannel[i] = -1;
1715 }// prevent division by 0
1717 if (portChannel[i] == -1) continue;
1719 Double_t dCOG = (Double_t)(rightCharge[i]-leftCharge[i])/centerCharge[i];
1721 Bool_t sign = (dCOG>=0.0) ? kFALSE : kTRUE;
1722 dCOG = (sign == kFALSE) ? dCOG : -dCOG; // AssignDouble doesn't allow for signed doubles
1723 dCOGAlu.AssignDouble(dCOG);
1724 Int_t iLUTpos = dCOGAlu.GetValue(); // steers position in LUT
1727 yrawAlu.AssignDouble(dCOG);
1728 Int_t iCOG = yrawAlu.GetValue();
1729 Int_t y = iCOG + fPosLUT[iLUTpos % 128]; // local position in pad-units
1730 yrawAlu.AssignFormatted(y); // 0<y<1
1731 yAlu = yrawAlu; // convert to 16 past-comma bits
1733 if(sign == kTRUE) yAlu.SetSign(-1); // buffer width of 9 bits; sign on real (not estimated) position
1734 xAlu.AssignInt(iT); // buffer width of 5 bits
1737 xxAlu = xAlu * xAlu; // buffer width of 10 bits -> fulfilled by x*x
1739 yyAlu = yAlu * yAlu; // buffer width of 16 bits
1741 xyAlu = xAlu * yAlu; // buffer width of 14 bits
1743 Int_t adc = portChannel[i]-1; // remapping! port-channel contains channel-nr. of inner adc's (1..19; mapped to 0..18)
1745 // calculate fit-sums recursively
1746 // interpretation of their bit-length is given as comment
1748 // be aware that the accuracy of the result of a calculation is always determined by the accuracy of the less accurate value
1750 XAlu.AssignFormatted(X[adc]);
1751 XAlu = XAlu + xAlu; // buffer width of 9 bits
1752 X[adc] = XAlu.GetValue();
1754 XXAlu.AssignFormatted(XX[adc]);
1755 XXAlu = XXAlu + xxAlu; // buffer width of 14 bits
1756 XX[adc] = XXAlu.GetValue();
1759 YAlu.AssignFormatted(-Y[adc]); // make sure that only positive values are assigned; sign-setting must be done by hand
1763 YAlu.AssignFormatted(Y[adc]);
1767 YAlu = YAlu + yAlu; // buffer width of 14 bits (8 past-comma);
1768 Y[adc] = YAlu.GetSignedValue();
1770 YYAlu.AssignFormatted(YY[adc]);
1771 YYAlu = YYAlu + yyAlu; // buffer width of 21 bits (16 past-comma)
1772 YY[adc] = YYAlu.GetValue();
1775 XYAlu.AssignFormatted(-XY[adc]);
1779 XYAlu.AssignFormatted(XY[adc]);
1783 XYAlu = XYAlu + xyAlu; // buffer allows 17 bits (8 past-comma)
1784 XY[adc] = XYAlu.GetSignedValue();
1786 N[adc] = N[adc] + 1;
1789 // accumulated charge
1790 qsumAlu.AssignFormatted(qsum[adc+1]); // qsum was not remapped!
1791 qtruncAlu = qsumAlu;
1793 if(iT>=tQS0 && iT<=tQE0){
1794 QT0Alu.AssignFormatted(QT0[adc]);
1795 QT0Alu = QT0Alu + qtruncAlu;
1796 QT0[adc] = QT0Alu.GetValue();
1797 //interpretation of QT0 as 12bit-value (all pre-comma); is this as it should be done?; buffer allows 15 Bit
1800 if(iT>=tQS1 && iT<=tQE1){
1801 QT1Alu.AssignFormatted(QT1[adc]);
1802 QT1Alu = QT1Alu + qtruncAlu;
1803 QT1[adc] = QT1Alu.GetValue();
1804 //interpretation of QT1 as 12bit-value; buffer allows 16 Bit
1808 // remapping is done!!
1814 // tracklet-assembly
1816 // put into AliTRDfeeParam and take care that values are in proper range
1817 const Int_t cTCL = 1; // left adc: number of hits; 8<=TCL<=31 (?? 1<=cTCL<+8 ??)
1818 const Int_t cTCT = 8; // joint number of hits; 8<=TCT<=31; note that according to TRAP manual this number cannot be lower than 8; however it should be adjustable to the number of hits in the fit time range (40%)
1820 Int_t mPair = 0; // marker for possible tracklet pairs
1821 Int_t* hitSum = new Int_t[fNADC-3];
1822 // hitSum[0] means: hit sum of remapped channels 0 and 1; hitSum[17]: 17 and 18;
1824 // check for all possible tracklet-pairs of adjacent channels (two are merged); mark the left channel of the chosen pairs
1825 for (Int_t iCol = 0; iCol < fNADC-3; iCol++) {
1826 hitSum[iCol] = N[iCol] + N[iCol+1];
1827 if ((N[iCol]>=cTCL) && (hitSum[iCol]>=cTCT)) {
1828 mPair |= 1; // mark as possible channel-pair
1835 List_t* selectPair = new List_t[fNADC-2]; // list with 18 elements (0..18) containing the left channel-nr and hit sums
1836 // selectPair[18] is starting list-element just for pointing
1837 for(Int_t k = 0; k<fNADC-2; k++){
1838 selectPair[k].next = NULL;
1839 selectPair[k].iadc = -1; // invalid adc
1840 selectPair[k].value = 0;
1847 // read marker and sort according to hit-sum
1849 Int_t adcL = 0; // left adc-channel-number (remapped)
1850 Int_t selNr = 0; // current number in list
1852 // insert marked channels into list and sort according to hit-sum
1853 while(adcL < fNADC-3 && selNr < fNADC-3){
1855 if( ((mPair>>((fNADC-4)-(adcL))) & 1) == 1) {
1856 selectPair[selNr].iadc = adcL;
1857 selectPair[selNr].value = hitSum[adcL];
1859 listLeft = &selectPair[fNADC-3];
1860 list = listLeft->next;
1863 while((list->next != NULL) && (selectPair[selNr].value <= list->value)){
1868 if(selectPair[selNr].value <= list->value){
1869 selectPair[selNr].next = list->next;
1870 list->next = &selectPair[selNr];
1873 listLeft->next = &selectPair[selNr];
1874 selectPair[selNr].next = list;
1879 listLeft->next = &selectPair[selNr];
1880 selectPair[selNr].next = list;
1888 //select up to 4 channels with maximum number of hits
1889 Int_t lpairChannel[4] = {-1,-1,-1,-1}; // save the left channel-numbers of pairs with most hit-sum
1890 Int_t rpairChannel[4] = {-1,-1,-1,-1}; // save the right channel, too; needed for detecting double tracklets
1891 list = &selectPair[fNADC-3];
1893 for (Int_t i = 0; i<4; i++) {
1894 if(list->next == NULL) continue;
1896 if(list->iadc == -1) continue;
1897 lpairChannel[i] = list->iadc; // channel number with selected hit
1898 rpairChannel[i] = lpairChannel[i]+1;
1901 // avoid submission of double tracklets
1902 for (Int_t i = 3; i>0; i--) {
1903 for (Int_t j = i-1; j>-1; j--) {
1904 if(lpairChannel[i] == rpairChannel[j]) {
1905 lpairChannel[i] = -1;
1906 rpairChannel[i] = -1;
1909 /* if(rpairChannel[i] == lpairChannel[j]) {
1910 lpairChannel[i] = -1;
1911 rpairChannel[i] = -1;
1917 // merging of the fit-sums of the remainig channels
1918 // assume same data-word-width as for fit-sums for 1 channel
1931 Int_t mMeanCharge[4];
1936 for (Int_t i = 0; i<4; i++){
1937 mADC[i] = -1; // set to invalid number
1951 oneAlu.AssignInt(1);
1952 one = oneAlu.GetValue(); // one with 8 past comma bits
1954 for (Int_t i = 0; i<4; i++){
1957 mADC[i] = lpairChannel[i]; // mapping of merged sums to left channel nr. (0,1->0; 1,2->1; ... 17,18->17)
1958 // the adc and pad-mapping should now be one to one: adc i is linked to pad i; TRAP-numbering
1959 Int_t madc = mADC[i];
1960 if (madc == -1) continue;
1962 YAlu.AssignInt(N[rpairChannel[i]]);
1963 Int_t wpad = YAlu.GetValue(); // enlarge hit counter of right channel by 8 past-comma bits; YAlu can have 5 pre-comma bits (values up to 63); hit counter<=nr of time bins (24)
1965 mN[i] = hitSum[madc];
1967 // don't merge fit sums in case of a stand-alone tracklet (consisting of only 1 channel); in that case only left channel makes up the fit sums
1968 if (N[madc+1] == 0) {
1969 mQT0[i] = QT0[madc];
1970 mQT1[i] = QT1[madc];
1975 // is it ok to do the size-checking for the merged fit-sums with the same format as for single-channel fit-sums?
1977 mQT0[i] = QT0[madc] + QT0[madc+1];
1978 QT0Alu.AssignFormatted(mQT0[i]);
1979 QT0Alu = QT0Alu; // size-check
1980 mQT0[i] = QT0Alu.GetValue(); // write back
1982 mQT1[i] = QT1[madc] + QT1[madc+1];
1983 QT1Alu.AssignFormatted(mQT1[i]);
1985 mQT1[i] = QT1Alu.GetValue();
1988 // calculate the mean charge in adc values; later to be replaced by electron likelihood
1989 mMeanCharge[i] = mQT0[i] + mQT1[i]; // total charge
1990 mMeanCharge[i] = mMeanCharge[i]>>2; // losing of accuracy; accounts for high mean charge
1991 // simulate LUT for 1/N; LUT is fed with the double-accurate pre-calculated value of 1/N; accuracy of entries has to be adjusted to real TRAP
1993 inverseNAlu.AssignDouble(invN);
1994 inverseN = inverseNAlu.GetValue();
1995 mMeanCharge[i] = mMeanCharge[i] * inverseN; // now to be interpreted with 8 past-comma bits
1996 TotalChargeAlu.AssignFormatted(mMeanCharge[i]);
1997 TotalChargeAlu = TotalChargeAlu;
1998 MeanChargeAlu = TotalChargeAlu;
1999 mMeanCharge[i] = MeanChargeAlu.GetValue();
2001 // this check is not necessary; it is just for efficiency reasons
2002 if (N[madc+1] == 0) {
2011 mX[i] = X[madc] + X[madc+1];
2012 XAlu.AssignFormatted(mX[i]);
2014 mX[i] = XAlu.GetValue();
2016 mXX[i] = XX[madc] + XX[madc+1];
2017 XXAlu.AssignFormatted(mXX[i]);
2019 mXX[i] = XXAlu.GetValue();
2022 mY[i] = Y[madc] + Y[madc+1] + wpad;
2024 YAlu.AssignFormatted(-mY[i]);
2028 YAlu.AssignFormatted(mY[i]);
2032 mY[i] = YAlu.GetSignedValue();
2034 mXY[i] = XY[madc] + XY[madc+1] + X[madc+1]*one; // multiplication by one to maintain the data format
2037 XYAlu.AssignFormatted(-mXY[i]);
2041 XYAlu.AssignFormatted(mXY[i]);
2045 mXY[i] = XYAlu.GetSignedValue();
2047 mYY[i] = YY[madc] + YY[madc+1] + 2*Y[madc+1]*one+ wpad*one;
2049 YYAlu.AssignFormatted(-mYY[i]);
2053 YYAlu.AssignFormatted(mYY[i]);
2058 mYY[i] = YYAlu.GetSignedValue();
2063 // calculation of offset and slope from the merged fit-sums;
2064 // YY is needed for some error measure only; still to be done
2065 // be aware that all values are relative values (scale: timebin-width; pad-width) and are integer values on special scale
2067 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2068 // !!important note: the offset is calculated from hits in the time bin range between tFS and tFE; it corresponds to the value at the height of the time bin tFS which does NOT need to correspond to the upper side of the drift !!
2069 // !!volume (cathode wire plane). The offset cannot be rescaled as long as it is unknown which is the first time bin that contains hits from the drift region and thus to which distance from the cathode plane tFS corresponds. !!
2070 // !!This has to be taken into account by the GTU. Furthermore a Lorentz correction might have to be applied to the offset (see below). !!
2071 // !!In this implementation it is assumed that no miscalibration containing changing drift velocities in the amplification region is used. !!
2072 // !!The corrections to the offset (e.g. no ExB correction applied as offset is supposed to be on top of drift region; however not at anode wire, so some inclination of drifting clusters due to Lorentz angle exists) are only !!
2073 // !!valid (in approximation) if tFS is close to the beginning of the drift region. !!
2074 // !!The slope however can be converted to a deflection length between electrode and cathode wire plane as it is clear that the drift region is sampled 20 times !!
2075 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2077 // which formats should be chosen?
2078 AliTRDtrapAlu denomAlu;
2079 denomAlu.Init(20,8);
2080 AliTRDtrapAlu numAlu;
2082 // is this enough pre-comma place? covers the range of the 13 bit-word of the transmitted offset
2083 // offset measured in coord. of left channel must be between -0.5 and 1.5; 14 pre-comma bits because numerator can be big
2085 for (Int_t i = 0; i<4; i++) {
2086 if (mADC[i] == -1) continue;
2088 Int_t num0 = (mN[i]*mXX[i]-mX[i]*mX[i]);
2090 denomAlu.AssignInt(-num0); // num0 does not have to be interpreted as having past-comma bits -> AssignInt
2091 denomAlu.SetSign(-1);
2094 denomAlu.AssignInt(num0);
2095 denomAlu.SetSign(1);
2098 Int_t num1 = mN[i]*mXY[i] - mX[i]*mY[i];
2100 numAlu.AssignFormatted(-num1); // value of num1 is already formatted to have 8 past-comma bits
2104 numAlu.AssignFormatted(num1);
2107 numAlu = numAlu/denomAlu;
2108 mSlope[i] = numAlu.GetSignedValue();
2110 Int_t num2 = mXX[i]*mY[i] - mX[i]*mXY[i];
2113 numAlu.AssignFormatted(-num2);
2117 numAlu.AssignFormatted(num2);
2121 numAlu = numAlu/denomAlu;
2124 mOffset[i] = numAlu.GetSignedValue();
2126 denomAlu.SetSign(1);
2129 //numAlu.AssignInt(mADC[i]+1); // according to TRAP-manual but trafo not to middle of chamber (0.5 channels away)
2130 numAlu.AssignDouble((Double_t)mADC[i] + 1.5); // numAlu has enough pre-comma place for that; correct trafo, best values
2131 mOffset[i] = mOffset[i] + numAlu.GetValue(); // transform offset to a coord.system relative to chip; +1 to avoid neg. values
2133 // up to here: adc-mapping according to TRAP manual and in line with pad-col mapping
2134 // reverse adc-counting to be again in line with the online mapping
2135 mADC[i] = fNADC - 4 - mADC[i]; // fNADC-4-mADC[i]: 0..17; remapping necessary;
2136 mADC[i] = mADC[i] + 2;
2137 // +2: mapping onto original ADC-online-counting: inner adc's corresponding to a chip's pasa: number 2..19
2140 // adc-counting is corresponding to online mapping; use AliTRDfeeParam::GetPadColFromADC to get the pad to which adc is connected;
2141 // pad-column mapping is reverse to adc-online mapping; TRAP adc-mapping is in line with pad-mapping (increase in same direction);
2143 // transform parameters to the local coordinate-system of a stack (used by GTU)
2144 AliTRDpadPlane* padPlane = fGeo->CreatePadPlane(fLayer,fStack);
2146 Double_t padWidthI = padPlane->GetWidthIPad()*10.0; // get values in cm; want them in mm
2147 //Double_t padWidthO = padPlane->GetWidthOPad()*10; // difference between outer pad-widths not included; in real TRAP??
2149 // difference between width of inner and outer pads of a row is not accounted for;
2151 Double_t magField = 0.4; // z-component of magnetic field in Tesla; adjust to current simulation!!; magnetic field can hardly be evaluated for the position of each mcm
2152 Double_t eCharge = 0.3; // unit charge in (GeV/c)/m*T
2153 Double_t ptMin = 2.3; // minimum transverse momentum (GeV/c); to be adjusted(?)
2155 Double_t granularityOffset = 0.160; // granularity for offset in mm
2156 Double_t granularitySlope = 0.140; // granularity for slope in mm
2158 // get the coordinates in SM-system; parameters:
2160 Double_t zPos = (padPlane->GetRowPos(fRow))*10.0; // z-position of the MCM; fRow is counted on a chamber; SM consists of 5
2161 // zPos is position of pad-borders;
2162 Double_t zOffset = 0.0;
2163 if ( fRow == 0 || fRow == 15 ) {
2164 zOffset = padPlane->GetLengthOPad();
2167 zOffset = padPlane->GetLengthIPad();
2169 zOffset = (-1.0) * zOffset/2.0;
2170 // turn zPos to be z-coordinate at middle of pad-row
2171 zPos = zPos + zOffset;
2174 Double_t xPos = 0.0; // x-position of the upper border of the drift-chamber of actual layer
2175 Int_t icol = 0; // column-number of adc-channel
2176 Double_t yPos[4]; // y-position of the pad to which ADC is connected
2177 Double_t dx = 30.0; // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry
2178 Double_t freqSample = fFeeParam->GetSamplingFrequency(); // retrieve the sampling frequency (10.019750 MHz)
2179 Double_t vdrift = fCal->GetVdriftAverage(fChaId); // averaged drift velocity for this detector (1.500000 cm/us)
2180 Int_t nrOfDriftTimeBins = Int_t(dx/10.0*freqSample/vdrift); // the number of time bins in the drift region (20)
2181 Int_t nrOfAmplTimeBins = 2; // the number of time bins between anode wire and cathode wires in ampl.region (3.5mm)(guess)(suppose v_drift+3.5cm/us there=>all clusters arrive at anode wire within one time bin (100ns))
2182 Int_t nrOfOffsetCorrTimeBins = tFS - nrOfAmplTimeBins - 1; // -1 is to be conservative; offset correction will not remove the shift but is supposed to improve it; if tFS = 5, 2 drift time bins before tFS are assumed
2183 if(nrOfOffsetCorrTimeBins < 0) nrOfOffsetCorrTimeBins = 0;// don't apply offset correction if no drift time bins before tFS can be assumed
2184 Double_t lorTan = fCal->GetOmegaTau(vdrift,magField); // tan of the Lorentz-angle for this detector; could be evaluated and set as a parameter for each mcm
2185 //Double_t lorAngle = 7.0; // Lorentz-angle in degrees
2186 Double_t tiltAngle = padPlane->GetTiltingAngle(); // sign-respecting tilting angle of pads in actual layer
2187 Double_t tiltTan = TMath::Tan(TMath::Pi()/180.0 * tiltAngle);
2188 //Double_t lorTan = TMath::Tan(TMath::Pi()/180.0 * lorAngle);
2190 Double_t alphaMax[4]; // maximum deflection from the direction to the primary vertex; granularity of hit pads
2191 Double_t slopeMin[4]; // local limits for the deflection
2192 Double_t slopeMax[4];
2193 Int_t mslopeMin[4]; // in granularity units; to be compared to mSlope[i]
2197 // x coord. of upper side of drift chambers in local SM-system (in mm)
2198 // obtained by evaluating the x-range of the hits; should be crosschecked; only drift, not amplification region taken into account (30mm);
2199 // the y-deflection is given as difference of y between lower and upper side of drift-chamber, not pad-plane;
2222 // calculation of offset-correction n:
2224 Int_t nCorrectOffset = (fRobPos % 2 == 0) ? ((fMcmPos % 4)) : ( 4 + (fMcmPos % 4));
2226 nCorrectOffset = (nCorrectOffset - 4)*18 - 1;
2227 if (nCorrectOffset < 0) {
2228 numAlu.AssignInt(-nCorrectOffset);
2232 numAlu.AssignInt(nCorrectOffset);
2235 nCorrectOffset = numAlu.GetSignedValue();
2237 // the Lorentz correction to the offset
2238 Double_t lorCorrectOffset = lorTan *(Double_t)nrOfOffsetCorrTimeBins*vdrift*10.0/freqSample; // Lorentz offset correction in mm
2241 lorCorrectOffset = lorCorrectOffset/padWidthI; // Lorentz correction in pad width units
2243 if(lorCorrectOffset < 0) {
2244 numAlu.AssignDouble(-lorCorrectOffset);
2248 numAlu.AssignDouble(lorCorrectOffset);
2252 Int_t mlorCorrectOffset = numAlu.GetSignedValue();
2255 Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0
2257 // calculation of slope-correction
2259 // this is only true for tracks coming (approx.) from primary vertex
2260 // everything is evaluated for a tracklet covering the whole drift chamber
2261 Double_t cCorrectSlope = (-lorTan*dx + zPos/xPos*dx*tiltTan)/granularitySlope;
2262 // Double_t cCorrectSlope = zPos/xPos*dx*tiltTan/granularitySlope;
2263 // zPos can be negative! for track from primary vertex: zOut-zIn > 0 <=> zPos > 0
2265 if (cCorrectSlope < 0) {
2266 numAlu.AssignDouble(-cCorrectSlope);
2270 numAlu.AssignDouble(cCorrectSlope);
2273 cCorrectSlope = numAlu.GetSignedValue();
2275 // convert slope to deflection between upper and lower drift-chamber position (slope is given in pad-unit/time-bins)
2276 // different pad-width of outer pads of a pad-plane not taken into account
2277 // note that the fit was only done in the range tFS to tFE, however this range does not need to cover the whole drift region (neither start nor end of it)
2278 // however the tracklets are supposed to be a fit in the drift region thus the linear function is stretched to fit the drift region of 30 mm
2281 Double_t mCorrectSlope = (Double_t)(nrOfDriftTimeBins)*padWidthI/granularitySlope; // >= 0.0
2283 AliTRDtrapAlu correctAlu;
2284 correctAlu.Init(20,8);
2286 AliTRDtrapAlu offsetAlu;
2287 offsetAlu.Init(13,0,-0x1000,0x0FFF); // 13 bit-word; 2-complement (1 sign-bit); asymmetric range
2289 AliTRDtrapAlu slopeAlu;
2290 slopeAlu.Init(7,0,-0x40,0x3F); // 7 bit-word; 2-complement (1 sign-bit);
2292 for (Int_t i = 0; i<4; i++) {
2294 if (mADC[i] == -1) continue;
2296 icol = fFeeParam->GetPadColFromADC(fRobPos,fMcmPos,mADC[i]); // be aware that mADC[i] contains the ADC-number according to online-mapping
2297 yPos[i] = (padPlane->GetColPos(icol))*10.0;
2302 correctAlu.AssignDouble(mCorrectOffset); // done because max. accuracy is 8 bit
2303 mCorrectOffset = correctAlu.GetValueWhole(); // cut offset correction to 8 past-comma bit accuracy
2304 mOffset[i] = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset - mlorCorrectOffset));
2305 //mOffset[i] = mOffset[i]*(-1); // adjust to direction of y-axes in online simulation
2307 if (mOffset[i] < 0) {
2308 numAlu.AssignFormatted(-mOffset[i]);
2312 numAlu.AssignFormatted(mOffset[i]);
2317 mOffset[i] = offsetAlu.GetSignedValue();
2322 correctAlu.AssignDouble(mCorrectSlope);
2323 mCorrectSlope = correctAlu.GetValueWhole();
2325 mSlope[i] = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope);
2327 if (mSlope[i] < 0) {
2328 numAlu.AssignFormatted(-mSlope[i]);
2332 numAlu.AssignFormatted(mSlope[i]);
2336 slopeAlu = numAlu; // here all past-comma values are cut, not rounded; alternatively add +0.5 before cutting (means rounding)
2337 mSlope[i] = slopeAlu.GetSignedValue();
2339 // local (LTU) limits for the deflection
2340 // ATan returns angles in radian
2341 alphaMax[i] = TMath::ASin(eCharge*magField/(2.0*ptMin)*(TMath::Sqrt(xPos*xPos + yPos[i]*yPos[i]))/1000.0); // /1000: mm->m
2342 slopeMin[i] = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) - alphaMax[i]))/granularitySlope;
2343 slopeMax[i] = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) + alphaMax[i]))/granularitySlope;
2345 if (slopeMin[i] < 0) {
2346 slopeAlu.AssignDouble(-slopeMin[i]);
2347 slopeAlu.SetSign(-1);
2350 slopeAlu.AssignDouble(slopeMin[i]);
2351 slopeAlu.SetSign(1);
2353 mslopeMin[i] = slopeAlu.GetSignedValue(); // the borders should lie inside the range of mSlope -> usage of slopeAlu again
2355 if (slopeMax[i] < 0) {
2356 slopeAlu.AssignDouble(-slopeMax[i]);
2357 slopeAlu.SetSign(-1);
2360 slopeAlu.AssignDouble(slopeMax[i]);
2361 slopeAlu.SetSign(1);
2363 mslopeMax[i] = slopeAlu.GetSignedValue();
2366 // suppress submission of tracks with low stiffness
2367 // put parameters in 32bit-word and submit (write to file as root-file; sort after SM, stack, layer, chamber)
2369 // sort tracklet-words in ascending y-order according to the offset (according to mADC would also be possible)
2370 // up to now they are sorted according to maximum hit sum
2371 // is the sorting really done in the TRAP-chip?
2373 Int_t order[4] = {-1,-1,-1,-1};
2374 Int_t wordnr = 0; // number of tracklet-words
2376 for(Int_t j = 0; j < fMaxTracklets; j++) {
2377 //if( mADC[j] == -1) continue;
2378 if( (mADC[j] == -1) || (mSlope[j] < mslopeMin[j]) || (mSlope[j] > mslopeMax[j])) continue; // this applies a pt-cut
2380 if( wordnr-1 == 0) {
2384 // wordnr-1>0, wordnr-1<4
2385 order[wordnr-1] = j;
2386 for( Int_t k = 0; k < wordnr-1; k++) {
2387 if( mOffset[j] < mOffset[order[k]] ) {
2388 for( Int_t l = wordnr-1; l > k; l-- ) {
2389 order[l] = order[l-1];
2398 // fill the bit-words in ascending order and without gaps
2399 UInt_t bitWord[4] = {0,0,0,0}; // attention: unsigned int to have real 32 bits (no 2-complement)
2400 for(Int_t j = 0; j < wordnr; j++) { // only "wordnr" tracklet-words
2401 //Bool_t rem1 = kTRUE;
2404 //bit-word is 2-complement and therefore without sign
2405 bitWord[j] = 1; // this is the starting 1 of the bit-word (at 33rd position); the 1 must be ignored
2413 /*printf("mean charge: %d\n",mMeanCharge[i]);
2414 printf("row: %d\n",fRow);
2415 printf("slope: %d\n",mSlope[i]);
2416 printf("pad position: %d\n",mOffset[i]);
2417 printf("channel: %d\n",mADC[i]);*/
2419 // electron probability (currently not implemented; the mean charge is just scaled)
2420 shift = (UInt_t)mMeanCharge[i];
2421 for(Int_t iBit = 0; iBit < 8; iBit++) {
2422 bitWord[j] = bitWord[j]<<1;
2423 bitWord[j] |= (shift>>(7-iBit))&1;
2428 shift = (UInt_t)fRow;
2429 for(Int_t iBit = 0; iBit < 4; iBit++) {
2430 bitWord[j] = bitWord[j]<<1;
2431 bitWord[j] |= (shift>>(3-iBit))&1;
2432 //printf("%d", (fRow>>(3-iBit))&1);
2435 // deflection length
2437 shift = (UInt_t)(-mSlope[i]);
2438 // shift2 is 2-complement of shift
2440 for(Int_t iBit = 1; iBit < 7; iBit++) {
2442 shift2 |= (1- (((shift)>>(6-iBit))&1) );
2443 //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2445 shift2 = shift2 + 1;
2447 for(Int_t iBit = 0; iBit < 7; iBit++) {
2448 bitWord[j] = bitWord[j]<<1;
2449 bitWord[j] |= (shift2>>(6-iBit))&1;
2450 //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2454 shift = (UInt_t)(mSlope[i]);
2455 bitWord[j] = bitWord[j]<<1;
2458 for(Int_t iBit = 1; iBit < 7; iBit++) {
2459 bitWord[j] = bitWord[j]<<1;
2460 bitWord[j] |= (shift>>(6-iBit))&1;
2461 //printf("%d",(mSlope[i]>>(6-iBit))&1);
2466 if(mOffset[i] < 0) {
2467 shift = (UInt_t)(-mOffset[i]);
2469 for(Int_t iBit = 1; iBit < 13; iBit++) {
2471 shift2 |= (1-(((shift)>>(12-iBit))&1));
2472 //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
2474 shift2 = shift2 + 1;
2476 for(Int_t iBit = 0; iBit < 13; iBit++) {
2477 bitWord[j] = bitWord[j]<<1;
2478 bitWord[j] |= (shift2>>(12-iBit))&1;
2479 //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2483 shift = (UInt_t)mOffset[i];
2484 bitWord[j] = bitWord[j]<<1;
2487 for(Int_t iBit = 1; iBit < 13; iBit++) {
2488 bitWord[j] = bitWord[j]<<1;
2489 bitWord[j] |= (shift>>(12-iBit))&1;
2490 //printf("%d",(mOffset[i]>>(12-iBit))&1);
2496 //printf("bitWord: %u\n",bitWord[j]);
2497 //printf("adc: %d\n",mADC[i]);
2498 fMCMT[j] = bitWord[j];
2517 delete [] selectPair;
2521 //if you want to activate the MC tracklet output, set fgkMCTrackletOutput=kTRUE in AliTRDfeeParam
2523 if (!fFeeParam->GetMCTrackletOutput())
2526 AliLog::SetClassDebugLevel("AliTRDmcmSim", 10);
2527 AliLog::SetFileOutput("../log/tracklet.log");
2529 // testing for wordnr in order to speed up the simulation
2533 UInt_t *trackletWord = new UInt_t[fMaxTracklets];
2534 Int_t *adcChannel = new Int_t[fMaxTracklets];
2535 Int_t *trackRef = new Int_t[fMaxTracklets];
2539 AliTRDdigitsManager *digman = new AliTRDdigitsManager();
2540 digman->ReadDigits(gAlice->GetRunLoader()->GetLoader("TRDLoader")->TreeD());
2541 digman->SetUseDictionaries(kTRUE);
2542 AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
2544 for (Int_t j = 0; j < fMaxTracklets; j++) {
2546 trackletWord[j] = 0;
2548 if (bitWord[j]!=0) {
2549 trackletWord[u] = bitWord[j];
2550 adcChannel[u] = mADC[i]; // mapping onto the original adc-array to be in line with the digits-adc-ordering (21 channels in total on 1 mcm, 18 belonging to pads); mADC[i] should be >-1 in case bitWord[i]>0
2552 // Finding label of MC track
2553 TH1F *hTrkRef = new TH1F("trackref", "trackref", 100000, 0, 100000);
2555 Int_t padcol = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u]);
2556 Int_t padcol_ngb = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u] - 1);
2557 Int_t padrow = 4 * (fRobPos / 2) + fMcmPos / 4;
2558 Int_t det = 30 * fSector + 6 * fStack + fLayer;
2559 for(Int_t iTimebin = feeParam->GetLinearFitStart(); iTimebin < feeParam->GetLinearFitEnd(); iTimebin++) {
2560 track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2561 track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2562 track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2563 hTrkRef->Fill(track[0]);
2564 if (track[1] != track[0] && track[1] != -1)
2565 hTrkRef->Fill(track[1]);
2566 if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2567 hTrkRef->Fill(track[2]);
2568 if (padcol_ngb >= 0) {
2569 track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2570 track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2571 track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2572 hTrkRef->Fill(track[0]);
2573 if (track[1] != track[0] && track[1] != -1)
2574 hTrkRef->Fill(track[1]);
2575 if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2576 hTrkRef->Fill(track[2]);
2579 trackRef[u] = hTrkRef->GetMaximumBin() - 1;
2585 AliDataLoader *dl = gAlice->GetRunLoader()->GetLoader("TRDLoader")->GetDataLoader("tracklets");
2587 AliError("Could not get the tracklets data loader!");
2590 TTree *trackletTree = dl->Tree();
2593 trackletTree = dl->Tree();
2595 AliTRDtrackletMCM *trkl = new AliTRDtrackletMCM();
2596 TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
2598 trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
2599 trkbranch->SetAddress(&trkl);
2601 for (Int_t iTracklet = 0; iTracklet < fMaxTracklets; iTracklet++) {
2602 if (trackletWord[iTracklet] == 0)
2604 trkl->SetTrackletWord(trackletWord[iTracklet]);
2605 trkl->SetDetector(30*fSector + 6*fStack + fLayer);
2606 trkl->SetROB(fRobPos);
2607 trkl->SetMCM(fMcmPos);
2608 trkl->SetLabel(trackRef[iTracklet]);
2609 trackletTree->Fill();
2612 dl->WriteData("OVERWRITE");
2615 delete [] trackletWord;
2616 delete [] adcChannel;
2621 // error measure for quality of fit (not necessarily needed for the trigger)
2622 // cluster quality threshold (not yet set)
2623 // electron probability
2625 //_____________________________________________________________________________________
2626 void AliTRDmcmSim::GeneratefZSM1Dim()
2629 // Generate the array fZSM1Dim necessary
2630 // for the method ProduceRawStream
2634 // Supressed zeros indicated by -1 in digits array
2635 for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ )
2637 for( Int_t it = 0 ; it < fNTimeBin ; it++ )
2640 if(fADCF[iadc][it]==-1) // If is a supressed value
2644 else // Not suppressed
2651 // Make the 1 dim projection
2652 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
2654 for( Int_t it = 0 ; it < fNTimeBin ; it++ )
2656 fZSM1Dim[iadc] &= fZSM[iadc][it];
2660 //_______________________________________________________________________________________
2661 void AliTRDmcmSim::CopyArrays()
2664 // Initialize filtered data array with raw data
2665 // Method added for internal consistency
2668 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
2670 for( Int_t it = 0 ; it < fNTimeBin ; it++ )
2672 fADCF[iadc][it] = fADCR[iadc][it];
2676 //_______________________________________________________________________________________
2677 void AliTRDmcmSim::StartfastZS(Int_t pads, Int_t timebins)
2680 // Initialize just the necessary elements to perform
2681 // the zero suppression in the digitizer
2684 fFeeParam = AliTRDfeeParam::Instance();
2685 fSimParam = AliTRDSimParam::Instance();
2687 fNTimeBin = timebins;
2691 fADCR = new Int_t *[fNADC];
2692 fADCF = new Int_t *[fNADC];
2693 fADCT = new Int_t *[fNADC];
2694 fZSM = new Int_t *[fNADC];
2695 fZSM1Dim = new Int_t [fNADC];
2696 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
2698 fADCR[iadc] = new Int_t[fNTimeBin];
2699 fADCF[iadc] = new Int_t[fNTimeBin];
2700 fADCT[iadc] = new Int_t[fNTimeBin];
2701 fZSM [iadc] = new Int_t[fNTimeBin];
2705 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
2707 for( Int_t it = 0 ; it < fNTimeBin ; it++ )
2709 fADCR[iadc][it] = 0;
2710 fADCF[iadc][it] = 0;
2711 fADCT[iadc][it] = -1;
2712 fZSM [iadc][it] = 1;
2717 fInitialized = kTRUE;
2719 //_______________________________________________________________________________________
2720 void AliTRDmcmSim::FlagDigitsArray(AliTRDarrayADC *tempdigs, Int_t valrow)
2723 // Modify the digits array to flag suppressed values
2726 for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ )
2728 for( Int_t it = 0 ; it < fNTimeBin ; it++ )
2730 if(fZSM[iadc][it]==1)
2732 tempdigs->SetData(valrow,iadc,it,-1);