1 /**************************************************************************
2 * Copyright(c) 1998-2000, 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 **************************************************************************/
19 Class for creating of the sumable digits and digits from MC data
21 The input : ideal signals (Hits->Diffusion->Attachment -Ideal signal)
22 The output: "raw digits"
25 1. Pad by pad gain map
27 3. The dead channels identified - zerro noise for corresponding pads
28 In this case the outpu equal zerro
37 #include <TObjArray.h>
39 #include <TDirectory.h>
40 #include <Riostream.h>
41 #include <TParameter.h>
43 #include "AliTPCDigitizer.h"
46 #include "AliTPCParam.h"
47 #include "AliTPCParamSR.h"
49 #include "AliLoader.h"
51 #include "AliDigitizationInput.h"
52 #include "AliSimDigits.h"
55 #include "AliTPCcalibDB.h"
56 #include "AliTPCCalPad.h"
57 #include "AliTPCCalROC.h"
58 #include "TTreeStream.h"
59 #include "AliTPCReconstructor.h"
60 #include <TGraphErrors.h>
65 ClassImp(AliTPCDigitizer)
67 //___________________________________________
68 AliTPCDigitizer::AliTPCDigitizer() :AliDigitizer(),fDebug(0), fDebugStreamer(0)
71 // Default ctor - don't use it
76 //___________________________________________
77 AliTPCDigitizer::AliTPCDigitizer(AliDigitizationInput* digInput)
78 :AliDigitizer(digInput),fDebug(0), fDebugStreamer(0)
81 // ctor which should be used
83 AliDebug(2,"(AliDigitizationInput* digInput) was processed");
84 if (AliTPCReconstructor::StreamLevel()>0) fDebugStreamer = new TTreeSRedirector("TPCDigitDebug.root","recreate");
88 //------------------------------------------------------------------------
89 AliTPCDigitizer::~AliTPCDigitizer()
92 if (fDebugStreamer) delete fDebugStreamer;
97 //------------------------------------------------------------------------
98 Bool_t AliTPCDigitizer::Init()
106 //------------------------------------------------------------------------
107 void AliTPCDigitizer::Digitize(Option_t* option)
109 //DigitizeFast(option);
110 DigitizeWithTailAndCrossTalk(option);
113 //------------------------------------------------------------------------
114 void AliTPCDigitizer::DigitizeFast(Option_t* option)
117 // merge input tree's with summable digits
118 //output stored in TreeTPCD
121 TString optionString = option;
122 if (!strcmp(optionString.Data(),"deb")) {
123 cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
126 //get detector and geometry
129 AliRunLoader *rl, *orl;
130 AliLoader *gime, *ogime;
134 Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
135 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
138 Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
144 AliTPC *pTPC = (AliTPC *) gAlice->GetModule("TPC");
145 AliTPCParam * param = pTPC->GetParam();
147 //sprintf(s,param->GetTitle());
148 snprintf(s,100,"%s",param->GetTitle());
149 //sprintf(ss,"75x40_100x60");
150 snprintf(ss,100,"75x40_100x60");
152 printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
154 param=new AliTPCParamSR();
157 //sprintf(ss,"75x40_100x60_150x60");
158 snprintf(ss,100,"75x40_100x60_150x60");
159 if(strcmp(s,ss)!=0) {
160 printf("No TPC parameters found...\n");
165 pTPC->GenerNoise(500000); //create table with noise
167 Int_t nInputs = fDigInput->GetNinputs();
168 Int_t * masks = new Int_t[nInputs];
169 for (Int_t i=0; i<nInputs;i++)
170 masks[i]= fDigInput->GetMask(i);
171 Short_t **pdig= new Short_t*[nInputs]; //pointers to the expanded digits array
172 Int_t **ptr= new Int_t*[nInputs]; //pointers to the expanded tracks array
173 Bool_t *active= new Bool_t[nInputs]; //flag for active input segments
176 //create digits array for given sectors
179 //create branch's in TPC treeD
180 orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
181 ogime = orl->GetLoader("TPCLoader");
182 TTree * tree = ogime->TreeD();
183 AliSimDigits * digrow = new AliSimDigits;
187 ogime->MakeTree("D");
188 tree = ogime->TreeD();
190 tree->Branch("Segment","AliSimDigits",&digrow);
192 AliSimDigits ** digarr = new AliSimDigits*[nInputs];
193 for (Int_t i1=0;i1<nInputs; i1++)
197 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
198 gime = rl->GetLoader("TPCLoader");
199 gime->LoadSDigits("read");
200 TTree * treear = gime->TreeS();
204 cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
205 <<" input "<< i1<<endl;
206 for (Int_t i2=0;i2<i1+1; i2++){
208 if(digarr[i2]) delete digarr[i2];
218 //sprintf(phname,"lhcphase%d",i1);
219 snprintf(phname,100,"lhcphase%d",i1);
220 TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
221 ->FindObject("lhcphase0");
223 cerr<<"AliTPCDigitizer: LHC phase not found in"
224 <<" input "<< i1<<endl;
225 for (Int_t i2=0;i2<i1+1; i2++){
226 if(digarr[i2]) delete digarr[i2];
235 tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
237 if (treear->GetIndex()==0)
238 treear->BuildIndex("fSegmentID","fSegmentID");
239 treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
247 param->SetZeroSup(2);
249 Int_t zerosup = param->GetZeroSup();
250 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
251 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
253 //Loop over segments of the TPC
255 for (Int_t segmentID=0; segmentID<param->GetNRowsTotal(); segmentID++)
257 Int_t sector, padRow;
258 if (!param->AdjustSectorRow(segmentID,sector,padRow))
260 cerr<<"AliTPC warning: invalid segment ID ! "<<segmentID<<endl;
263 AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector); // pad gains per given sector
264 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector); // noise per given sector
265 digrow->SetID(segmentID);
270 Bool_t digitize = kFALSE;
271 for (Int_t i=0;i<nInputs; i++)
274 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
275 gime = rl->GetLoader("TPCLoader");
277 if (gime->TreeS()->GetEntryWithIndex(segmentID,segmentID) >= 0) {
278 digarr[i]->ExpandBuffer();
279 digarr[i]->ExpandTrackBuffer();
280 nTimeBins = digarr[i]->GetNRows();
281 nPads = digarr[i]->GetNCols();
283 if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
287 if (GetRegionOfInterest() && !digitize) break;
289 if (!digitize) continue;
291 digrow->Allocate(nTimeBins,nPads);
292 digrow->AllocateTrack(3);
295 Int_t label[1000]; //stack for 300 events
298 Int_t nElems = nTimeBins*nPads;
300 for (Int_t i=0;i<nInputs; i++)
302 pdig[i] = digarr[i]->GetDigits();
303 ptr[i] = digarr[i]->GetTracks();
306 Short_t *pdig1= digrow->GetDigits();
307 Int_t *ptr1= digrow->GetTracks() ;
311 for (Int_t elem=0;elem<nElems; elem++)
317 for (Int_t i=0;i<nInputs; i++) if (active[i])
319 // q += digarr[i]->GetDigitFast(rows,col);
322 for (Int_t tr=0;tr<3;tr++)
324 // Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
325 Int_t lab = ptr[i][tr*nElems];
326 if ( (lab > 1) && *(pdig[i])>zerosup)
328 label[labptr]=lab+masks[i];
335 q/=16.; //conversion factor
336 Float_t gain = gainROC->GetValue(padRow,elem/nTimeBins); // get gain for given - pad-row pad
338 //printf("problem\n");
341 Float_t noisePad = noiseROC->GetValue(padRow,elem/nTimeBins);
342 // Float_t noise = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
343 Float_t noise = pTPC->GetNoise();
348 if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
349 //digrow->SetDigitFast((Short_t)q,rows,col);
351 for (Int_t tr=0;tr<3;tr++)
354 ptr1[tr*nElems] = label[tr];
363 digrow->GlitchFilter();
365 digrow->CompresBuffer(1,zerosup);
366 digrow->CompresTrackBuffer(1);
368 if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
369 } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
372 orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
373 ogime = orl->GetLoader("TPCLoader");
374 ogime->WriteDigits("OVERWRITE");
376 //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
379 for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
389 //------------------------------------------------------------------------
390 void AliTPCDigitizer::DigitizeSave(Option_t* option)
393 // Merge input tree's with summable digits
394 // Output digits stored in TreeTPCD
396 // Not active for long time.
397 // Before adding modification (for ion tail calucation and for the crorsstalk) it should be
398 // checked one by one with currenlty used AliTPCDigitizer::DigitizeFast
400 TString optionString = option;
401 if (!strcmp(optionString.Data(),"deb")) {
402 cout<<"AliTPCDigitizer::Digitize: called with option deb "<<endl;
405 //get detector and geometry
406 AliRunLoader *rl, *orl;
407 AliLoader *gime, *ogime;
410 orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
411 ogime = orl->GetLoader("TPCLoader");
413 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
414 //gime = rl->GetLoader("TPCLoader");
415 rl->GetLoader("TPCLoader");
417 AliRun* alirun = rl->GetAliRun();
419 AliTPC *pTPC = (AliTPC *) alirun->GetModule("TPC");
420 AliTPCParam * param = pTPC->GetParam();
421 pTPC->GenerNoise(500000); //create teble with noise
422 printf("noise %f \n", param->GetNoise()*param->GetNoiseNormFac());
424 Int_t nInputs = fDigInput->GetNinputs();
425 // stupid protection...
426 if (nInputs <= 0) return;
428 Int_t * masks = new Int_t[nInputs];
429 for (Int_t i=0; i<nInputs;i++)
430 masks[i]= fDigInput->GetMask(i);
432 AliSimDigits ** digarr = new AliSimDigits*[nInputs];
433 for(Int_t ii=0;ii<nInputs;ii++) digarr[ii]=0;
435 for (Int_t i1=0;i1<nInputs; i1++)
439 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
440 gime = rl->GetLoader("TPCLoader");
442 TTree * treear = gime->TreeS();
445 cerr<<" TPC - not existing input = \n"<<i1<<" ";
447 for(Int_t i=0; i<nInputs; i++) delete digarr[i];
452 TBranch * br = treear->GetBranch("fSegmentID");
453 if (br) br->GetFile()->cd();
454 treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
457 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
458 gime = rl->GetLoader("TPCLoader");
459 Stat_t nentries = gime->TreeS()->GetEntries();
462 //create branch's in TPC treeD
463 AliSimDigits * digrow = new AliSimDigits;
464 TTree * tree = ogime->TreeD();
466 tree->Branch("Segment","AliSimDigits",&digrow);
467 param->SetZeroSup(2);
469 Int_t zerosup = param->GetZeroSup();
470 //Loop over segments of the TPC
472 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
473 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
474 for (Int_t n=0; n<nentries; n++) {
475 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
476 gime = rl->GetLoader("TPCLoader");
477 gime->TreeS()->GetEvent(n);
479 digarr[0]->ExpandBuffer();
480 digarr[0]->ExpandTrackBuffer();
483 for (Int_t i=1;i<nInputs; i++){
484 // fDigInput->GetInputTreeTPCS(i)->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
485 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
486 gime = rl->GetLoader("TPCLoader");
487 gime->TreeS()->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
488 digarr[i]->ExpandBuffer();
489 digarr[i]->ExpandTrackBuffer();
490 if ((digarr[0]->GetID()-digarr[i]->GetID())>0)
495 Int_t sector, padRow;
496 if (!param->AdjustSectorRow(digarr[0]->GetID(),sector,padRow)) {
497 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
501 AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector); // pad gains per given sector
502 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector); // noise per given sector
503 digrow->SetID(digarr[0]->GetID());
505 Int_t nTimeBins = digarr[0]->GetNRows();
506 Int_t nPads = digarr[0]->GetNCols();
507 digrow->Allocate(nTimeBins,nPads);
508 digrow->AllocateTrack(3);
511 Int_t label[1000]; //stack for 300 events
516 for (Int_t iTimeBin=0;iTimeBin<nTimeBins; iTimeBin++){ // iTimeBin
517 for (Int_t iPad=0;iPad<nPads; iPad++){ // pad
522 for (Int_t i=0;i<nInputs; i++){
523 q += digarr[i]->GetDigitFast(iTimeBin,iPad);
526 for (Int_t tr=0;tr<3;tr++) {
527 Int_t lab = digarr[i]->GetTrackIDFast(iTimeBin,iPad,tr);
528 //Int_t lab = ptr[i][tr*nElems];
530 label[labptr]=lab+masks[i];
538 q/=16.; //conversion factor
539 // Float_t noise = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
540 Float_t gain = gainROC->GetValue(padRow,iPad);
542 Float_t noisePad = noiseROC->GetValue(padRow, iPad);
544 Float_t noise = pTPC->GetNoise();
547 // here we can get digits from past and add signal
550 //for (Int_t jTimeBin=0; jTimeBin<iTimeBin; jTimeBin++)
557 if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
558 digrow->SetDigitFast((Short_t)q,iTimeBin,iPad);
559 // *pdig1 =Short_t(q);
560 for (Int_t tr=0;tr<3;tr++){
562 ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],iTimeBin,iPad,tr);
563 //ptr1[tr*nElems] = label[tr];
565 // ((AliSimDigits*)digrow)->SetTrackIDFast(-1,iTimeBin,iPad,tr);
566 // ptr1[tr*nElems] = 1;
574 digrow->CompresBuffer(1,zerosup);
575 digrow->CompresTrackBuffer(1);
577 if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
579 // printf("end TPC merging - end -Tree %s\t%p\n",fDigInput->GetInputTreeH(0)->GetName(),fDigInput->GetInputTreeH(0)->GetListOfBranches()->At(3));
580 //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
581 ogime->WriteDigits("OVERWRITE");
583 for (Int_t i=1;i<nInputs; i++)
585 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
586 gime = rl->GetLoader("TPCLoader");
587 gime->UnloadSDigits();
589 ogime->UnloadDigits();
592 for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
603 //------------------------------------------------------------------------
604 void AliTPCDigitizer::DigitizeWithTailAndCrossTalk(Option_t* option)
606 // Modified version of the digitization function
607 // Modification: adding the ion tail and crosstalk:
609 // pcstream used in order to visually inspect data
612 // Crosstalk simulation:
613 // 1.) Calculate per time bin mean charge (per pad) within anode wire segment
614 // 2.) Subsract for the clusters at given time bin fraction of (mean charge) normalized by add hoc constant
615 // AliTPCRecoParam::GetCrosstalkCorrection() (0 if not crosstalk, 1 if ideal crosstalk)
616 // for simplicity we are assuming that wire segents are related to pad-rows
617 // Wire segmentationn is obtatined from the
618 // AliTPCParam::GetWireSegment(Int_t sector, Int_t row); // to be implemented
619 // AliTPCParam::GetNPadsPerSegment(Int_t segmentID); // to be implemented
620 // 3.) Substract form the signal contribution from the previous tracks - Ion tail in case speified in the AliTPCRecoParam
621 // AliTPCRecoParam::GetUseIonTailCorrection()
623 // Ion tail simulation:
624 // 1.) Needs signal from pad+-1, taking signal from history
625 // merge input tree's with summable digits
626 // output stored in TreeTPCD
628 AliTPCcalibDB* const calib=AliTPCcalibDB::Instance();
629 AliTPCRecoParam *recoParam = calib->GetRecoParam(0);
630 AliDebug(1, Form(" recoParam->GetCrosstalkCorrection() = %f", recoParam->GetCrosstalkCorrection()));
631 AliDebug(1,Form(" recoParam->GetUseIonTailCorrection() = %d ", recoParam->GetUseIonTailCorrection()));
635 TString optionString = option;
636 if (!strcmp(optionString.Data(),"deb")) {
637 cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
641 // ======== get detector and geometry =======
643 AliRunLoader *rl, *orl;
644 AliLoader *gime, *ogime;
648 Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
649 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
652 Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
658 AliTPC *pTPC = (AliTPC *) gAlice->GetModule("TPC");
659 AliTPCParam * param = pTPC->GetParam();
661 //sprintf(s,param->GetTitle());
662 snprintf(s,100,"%s",param->GetTitle());
663 //sprintf(ss,"75x40_100x60");
664 snprintf(ss,100,"75x40_100x60");
666 printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
668 param=new AliTPCParamSR();
670 //sprintf(ss,"75x40_100x60_150x60");
671 snprintf(ss,100,"75x40_100x60_150x60");
672 if(strcmp(s,ss)!=0) {
673 printf("No TPC parameters found...\n");
678 pTPC->GenerNoise(500000); //create table with noise
680 Int_t nInputs = fDigInput->GetNinputs();
681 Int_t * masks = new Int_t[nInputs];
682 for (Int_t i=0; i<nInputs;i++)
683 masks[i]= fDigInput->GetMask(i);
684 Short_t **pdig= new Short_t*[nInputs]; //pointers to the expanded digits array
685 Int_t **ptr= new Int_t*[nInputs]; //pointers to the expanded tracks array
686 Bool_t *active= new Bool_t[nInputs]; //flag for active input segments
689 //create digits array for given sectors
692 //create branch's in TPC treeD
693 orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
694 ogime = orl->GetLoader("TPCLoader");
695 TTree * tree = ogime->TreeD();
696 AliSimDigits * digrow = new AliSimDigits;
700 ogime->MakeTree("D");
701 tree = ogime->TreeD();
703 tree->Branch("Segment","AliSimDigits",&digrow);
705 AliSimDigits ** digarr = new AliSimDigits*[nInputs];
706 for (Int_t i1=0;i1<nInputs; i1++)
710 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
711 gime = rl->GetLoader("TPCLoader");
712 gime->LoadSDigits("read");
713 TTree * treear = gime->TreeS();
717 cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
718 <<" input "<< i1<<endl;
719 for (Int_t i2=0;i2<i1+1; i2++){
721 if(digarr[i2]) delete digarr[i2];
731 //sprintf(phname,"lhcphase%d",i1);
732 snprintf(phname,100,"lhcphase%d",i1);
733 TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
734 ->FindObject("lhcphase0");
737 cerr<<"AliTPCDigitizer: LHC phase not found in"
738 <<" input "<< i1<<endl;
739 for (Int_t i2=0;i2<i1+1; i2++){
740 if(digarr[i2]) delete digarr[i2];
749 tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
751 if (treear->GetIndex()==0)
752 treear->BuildIndex("fSegmentID","fSegmentID");
753 treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
760 // zero supp, take gain and noise map of TPC from OCDB
761 param->SetZeroSup(2);
762 Int_t zerosup = param->GetZeroSup();
763 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
764 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
768 // Cache the ion tail objects form OCDB
771 TObjArray *ionTailArr = (TObjArray*)AliTPCcalibDB::Instance()->GetIonTailArray();
772 if (!ionTailArr) {AliFatal("TPC - Missing IonTail OCDB object");}
773 // TObject *rocFactorIROC = ionTailArr->FindObject("factorIROC");
774 // TObject *rocFactorOROC = ionTailArr->FindObject("factorOROC");
775 // Float_t factorIROC = (atof(rocFactorIROC->GetTitle()));
776 // Float_t factorOROC = (atof(rocFactorOROC->GetTitle()));
777 Int_t nIonTailBins =0;
778 TObjArray timeResFunc(nROCs);
779 for (Int_t isec = 0;isec<nROCs;isec++){ //loop overs sectors
780 // Array of TGraphErrors for a given sector
781 TGraphErrors ** graphRes = new TGraphErrors *[20];
782 Float_t * trfIndexArr = new Float_t[20];
783 for (Int_t icache=0; icache<20; icache++)
785 graphRes[icache] = NULL;
786 trfIndexArr[icache] = 0;
788 if (!AliTPCcalibDB::Instance()->GetTailcancelationGraphs(isec,graphRes,trfIndexArr)) continue;
790 // fill all TGraphErrors of trfs (time response function) of a given sector to a TObjArray
791 TObjArray *timeResArr = new TObjArray(20); timeResArr -> SetOwner(kTRUE);
792 for (Int_t ires = 0;ires<20;ires++) timeResArr->AddAt(graphRes[ires],ires);
793 timeResFunc.AddAt(timeResArr,isec); // Fill all trfs into a single TObjArray
794 nIonTailBins = graphRes[3]->GetN();
798 // 1.) Make first loop to calculate mean amplitude per pad per segment for cross talk
801 TObjArray crossTalkSignalArray(nROCs); // for 36 sectors
802 TVectorD * qTotSectorOld = new TVectorD(nROCs);
803 Float_t qTotTPC = 0.;
804 Float_t qTotPerSector = 0.;
805 Int_t nTimeBinsAll = 1100;
806 Int_t nWireSegments=11;
807 // 1.a) crorstalk matrix initialization
808 for (Int_t sector=0; sector<nROCs; sector++){
809 TMatrixD *pcrossTalkSignal = new TMatrixD(nWireSegments,nTimeBinsAll);
810 for (Int_t imatrix = 0; imatrix<11; imatrix++)
811 for (Int_t jmatrix = 0; jmatrix<nTimeBinsAll; jmatrix++){
812 (*pcrossTalkSignal)[imatrix][jmatrix]=0.;
814 crossTalkSignalArray.AddAt(pcrossTalkSignal,sector);
817 // main loop over rows of whole TPC
818 for (Int_t globalRowID=0; globalRowID<param->GetNRowsTotal(); globalRowID++) {
819 Int_t sector, padRow;
820 if (!param->AdjustSectorRow(globalRowID,sector,padRow)) {
821 cerr<<"AliTPC warning: invalid segment ID ! "<<globalRowID<<endl;
824 // Calculate number of pads in a anode wire segment for normalization
825 Int_t wireSegmentID = param->GetWireSegment(sector,padRow);
826 Float_t nPadsPerSegment = (Float_t)(param->GetNPadsPerSegment(wireSegmentID));
827 // structure with mean signal per pad to be filled for each timebin in first loop (11 anodeWireSegment and 1100 timebin)
828 TMatrixD &crossTalkSignal = *((TMatrixD*)crossTalkSignalArray.At(sector));
829 AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector); // pad gains per given sector
830 // digrow->SetID(globalRowID);
833 Bool_t digitize = kFALSE;
834 for (Int_t i=0;i<nInputs; i++){ //here we can have more than one input - merging of separate events, signal1+signal2+background
835 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
836 gime = rl->GetLoader("TPCLoader");
837 if (gime->TreeS()->GetEntryWithIndex(globalRowID,globalRowID) >= 0) {
838 digarr[i]->ExpandBuffer();
839 digarr[i]->ExpandTrackBuffer();
840 nTimeBins = digarr[i]->GetNRows();
841 nPads = digarr[i]->GetNCols();
843 if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
847 if (GetRegionOfInterest() && !digitize) break;
849 if (!digitize) continue;
850 //digrow->Allocate(nTimeBins,nPads);
853 Int_t nElems = nTimeBins*nPads; // element is a unit of a given row's pad-timebin space
854 for (Int_t i=0;i<nInputs; i++)
856 pdig[i] = digarr[i]->GetDigits();
859 // loop over elements i.e pad-timebin space of a row, "padNumber=elem/nTimeBins", "timeBin=elem%nTimeBins"
861 for (Int_t elem=0;elem<nElems; elem++) {
865 for (Int_t i=0;i<nInputs; i++) if (active[i]) {
870 Int_t padNumber = elem/nTimeBins;
871 Int_t timeBin = elem%nTimeBins;
872 Float_t gain = gainROC->GetValue(padRow,padNumber); // get gain for given - pad-row pad
874 crossTalkSignal[wireSegmentID][timeBin]+= q/nPadsPerSegment; // Qtot per segment for a given timebin
875 qTotSectorOld -> GetMatrixArray()[sector] += q; // Qtot for each sector
876 qTotTPC += q; // Qtot for whole TPC
878 } // end of global row loop
881 // 1.b) Dump the content of the crossTalk signal to the debug stremer - to be corealted later with the same crosstalk correction
882 // assumed during reconstruction
884 if (AliTPCReconstructor::StreamLevel()==1) {
886 // dump the crosstalk matrices to tree for further investigation
887 // a.) to estimate fluctuation of pedestal in indiviula wire segments
888 // b.) to check correlation between regions
889 // c.) to check relative conribution of signal below threshold to crosstalk
890 for (Int_t isector=0; isector<nROCs; isector++){ //set all ellemts of crosstalk matrix to 0
891 TMatrixD * crossTalkMatrix = (TMatrixD*)crossTalkSignalArray.At(isector);
892 //TMatrixD * crossTalkMatrixBelow = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs);
893 TVectorD vecAll(crossTalkMatrix->GetNrows());
894 //TVectorD vecBelow(crossTalkMatrix->GetNrows());
896 for (Int_t itime=0; itime<crossTalkMatrix->GetNcols(); itime++){
897 for (Int_t iwire=0; iwire<crossTalkMatrix->GetNrows(); iwire++){
898 vecAll[iwire]=(*crossTalkMatrix)(iwire,itime);
899 //vecBelow[iwire]=(*crossTalkMatrixBelow)(iwire,itime);
901 (*fDebugStreamer)<<"crosstalkMatrix"<<
904 "vecAll.="<<&vecAll<< // crosstalk charge + charge below threshold
905 //"vecBelow.="<<&vecBelow<< // crosstalk contribution from signal below threshold
915 // 2.) Loop over segments (padrows) of the TPC
917 for (Int_t globalRowID=0; globalRowID<param->GetNRowsTotal(); globalRowID++) {
918 Int_t sector, padRow;
919 if (!param->AdjustSectorRow(globalRowID,sector,padRow)) {
920 cerr<<"AliTPC warning: invalid segment ID ! "<<globalRowID<<endl;
923 TObjArray *arrTRF = (TObjArray*)timeResFunc.At(sector);
924 // TGraphErrors *graphTRF = (TGraphErrors*)arrTRF->At(1);
925 Int_t wireSegmentID = param->GetWireSegment(sector,padRow);
926 Float_t nPadsPerSegment = (Float_t)(param->GetNPadsPerSegment(wireSegmentID));
927 // const Float_t ampfactor = (sector<36)?factorIROC:factorOROC; // factor for the iontail which is ROC type dependent
928 AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector); // pad gains per given sector
929 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector); // noise per given sector
930 digrow->SetID(globalRowID);
933 Bool_t digitize = kFALSE;
934 for (Int_t i=0;i<nInputs; i++) { //here we can have more than one input - merging of separate events, signal1+signal2+background
935 rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
936 gime = rl->GetLoader("TPCLoader");
937 if (gime->TreeS()->GetEntryWithIndex(globalRowID,globalRowID) >= 0) {
938 digarr[i]->ExpandBuffer();
939 digarr[i]->ExpandTrackBuffer();
940 nTimeBins = digarr[i]->GetNRows();
941 nPads = digarr[i]->GetNCols();
943 if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
947 if (GetRegionOfInterest() && !digitize) break;
949 if (!digitize) continue;
951 digrow->Allocate(nTimeBins,nPads);
952 digrow->AllocateTrack(3);
957 Float_t qIonTail = 0.;
959 Int_t label[1000]; //stack for 300 events
961 Int_t nElems = nTimeBins*nPads; // element is a unit of a given row's pad-timebin space
962 for (Int_t i=0;i<nInputs; i++)
964 pdig[i] = digarr[i]->GetDigits();
965 ptr[i] = digarr[i]->GetTracks();
967 Short_t *pdig1= digrow->GetDigits();
968 Int_t *ptr1= digrow->GetTracks() ;
969 // loop over elements i.e pad-timebin space of a row
970 for (Int_t elem=0;elem<nElems; elem++) {
974 for (Int_t i=0;i<nInputs; i++) if (active[i]){
976 for (Int_t tr=0;tr<3;tr++) {
977 Int_t lab = ptr[i][tr*nElems];
978 if ( (lab > 1) && *(pdig[i])>zerosup) {
979 label[labptr]=lab+masks[i];
986 Int_t padNumber = elem/nTimeBins;
987 Int_t timeBin = elem%nTimeBins;
988 localPad = padNumber-nPads/2;
990 Float_t gain = gainROC->GetValue(padRow,padNumber); // get gain for given - pad-row pad
992 //printf("problem\n");
996 Float_t noisePad = noiseROC->GetValue(padRow,padNumber);
997 Float_t noise = pTPC->GetNoise()*noisePad;
998 if ( (q/16.+noise)> zerosup){
999 // Crosstalk correction
1000 qXtalk = (*(TMatrixD*)crossTalkSignalArray.At(sector))[wireSegmentID][timeBin];
1001 qTotPerSector = qTotSectorOld -> GetMatrixArray()[sector];
1003 // Ion tail correction: being elem=padNumber*nTimeBins+timeBin;
1004 Int_t lowerElem=elem-nIonTailBins;
1005 Int_t zeroElem =(elem/nTimeBins)*nTimeBins;
1006 if (zeroElem<0) zeroElem=0;
1007 if (lowerElem<zeroElem) lowerElem=zeroElem;
1010 if (q>0 && recoParam->GetUseIonTailCorrection()){
1011 for (Int_t i=0;i<nInputs; i++) if (active[i]){
1012 Short_t *pdigC= digarr[i]->GetDigits();
1013 if (padNumber==0) continue;
1014 if (padNumber>=nPads-1) continue;
1015 for (Int_t dpad=-1; dpad<=1; dpad++){ // calculate iontail due signals from neigborhood pads
1016 for (Int_t celem=elem-1; celem>lowerElem; celem--){
1017 Int_t celemPad=celem+dpad*nTimeBins;
1018 Double_t qCElem=pdigC[celemPad];
1019 if ( qCElem<=0) continue;
1020 //here we substract ion tail
1022 if (celemPad-nTimeBins>nTimeBins && celemPad+nTimeBins<nElems){ // COG calculation in respect to current pad in pad units
1023 Double_t sumAmp=pdigC[celemPad-nTimeBins]+pdigC[celemPad]+pdigC[celemPad+nTimeBins];
1024 COG=(-1.0*pdigC[celemPad-nTimeBins]+pdigC[celemPad+nTimeBins])/sumAmp;
1026 Int_t indexTRFPRF = (TMath::Nint(TMath::Abs(COG*10.))%20);
1027 TGraphErrors *graphTRFPRF = (TGraphErrors*)arrTRF->At(indexTRFPRF);
1028 if (graphTRFPRF==NULL) continue;
1029 // here we should get index and point of TRF corresponding to given COG position
1030 if (graphTRFPRF->GetY()[elem-celem]<0)qIonTail+=qCElem*graphTRFPRF->GetY()[elem-celem];
1037 q -= qXtalk*recoParam->GetCrosstalkCorrection();
1039 q/=16.; //conversion factor
1041 q=TMath::Nint(q); // round to the nearest integer
1044 // fill info for degugging
1045 if (AliTPCReconstructor::StreamLevel()==1 && qOrig > zerosup ) {
1046 TTreeSRedirector &cstream = *fDebugStreamer;
1047 UInt_t uid = AliTPCROC::GetTPCUniqueID(sector, padRow, padNumber);
1048 cstream <<"ionTailXtalk"<<
1049 "uid="<<uid<< // globla unique identifier
1050 "sector="<< sector<<
1051 "globalRowID="<<globalRowID<<
1052 "padRow="<< padRow<< //pad row
1053 "wireSegmentID="<< wireSegmentID<< //wire segment 0-11, 0-3 in IROC 4-10 in OROC
1054 "localPad="<<localPad<< // pad number -npads/2 .. npads/2
1055 "padNumber="<<padNumber<< // pad number 0..npads
1056 "timeBin="<< timeBin<< // time bin
1057 "nPadsPerSegment="<<nPadsPerSegment<<// number of pads per wire segment
1058 "qTotPerSector="<<qTotPerSector<< // total charge in sector
1060 "qTotTPC="<<qTotTPC<< // acumulated charge without crosstalk and ion tail in full TPC
1061 "qOrig="<< qOrig<< // charge in given pad-row,pad,time-bin
1062 "q="<<q<< // q=qOrig-qXtalk-qIonTail - to check sign of the effects
1063 "qXtalk="<<qXtalk<< // crosstal contribtion at given position
1064 "qIonTail="<<qIonTail<< // ion tail cotribution from past signal
1066 }// dump the results to the debug streamer if in debug mode
1069 if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
1070 //digrow->SetDigitFast((Short_t)q,rows,col);
1072 for (Int_t tr=0;tr<3;tr++)
1075 ptr1[tr*nElems] = label[tr];
1082 if (AliTPCReconstructor::StreamLevel()==1 && qOrig > zerosup ) {
1083 cout << " sector = " << sector << " row = " << padRow << " localPad = " << localPad << " SegmentID = " << wireSegmentID << " nPadsPerSegment = " << nPadsPerSegment << endl;
1084 cout << " qXtalk = " << qXtalk ;
1085 cout << " qOrig = " << qOrig ;
1086 cout << " q = " << q ;
1087 cout << " qsec = " << qTotPerSector ;
1088 cout << " qTotTPC "<< qTotTPC << endl;
1093 digrow->GlitchFilter();
1095 digrow->CompresBuffer(1,zerosup);
1096 digrow->CompresTrackBuffer(1);
1098 if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
1099 } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
1102 orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
1103 ogime = orl->GetLoader("TPCLoader");
1104 ogime->WriteDigits("OVERWRITE");
1106 //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
1109 for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];