]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/Rec/AliTPCclustererKr.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
end-of-line normalization
[u/mrichter/AliRoot.git] / TPC / Rec / AliTPCclustererKr.cxx
CommitLineData
a65a7e70 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id: AliTPCclustererKr.cxx,v 1.7 2008/02/06 17:24:53 matyja Exp $ */
17
18//-----------------------------------------------------------------
19// Implementation of the TPC Kr cluster class
20//
21// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl
22//-----------------------------------------------------------------
23
24/*
25Instruction - how to use that:
26There are two macros prepared. One is for preparing clusters from MC
27samples: FindKrClusters.C. The output is kept in TPC.RecPoints.root.
28The other macro is prepared for data analysis: FindKrClustersRaw.C.
29The output is created for each processed file in root file named adc.root.
30For each data subsample the same named file is created. So be careful
31do not overwrite them.
32
33Additional selection criteria to select the GOLD cluster
34Example:
35// open file with clusters
36TFile f("Krypton.root");
37TTree * tree = (TTree*)f.Get("Kr")
38TCut cutR0("cutR0","fADCcluster/fSize<100"); // adjust it according v seetings -
39TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal
40TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.2"); // digital noise removal
41TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal
42TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks
43TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1;
44This values are typical values to be applied in selectors
45
46
47*
48**** MC ****
49*
50
51To run clusterizaton for MC type:
52.x FindKrClusters.C
53
54If you don't want to use the standard selection criteria then you
55have to do following:
56
57// load the standard setup
58AliRunLoader* rl = AliRunLoader::Open("galice.root");
59AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
60tpcl->LoadDigits();
61rl->LoadgAlice();
62gAlice=rl->GetAliRun();
63TDirectory *cwd = gDirectory;
64AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC");
65Int_t ver = tpc->IsVersion();
66rl->CdGAFile();
67AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");
68AliTPCDigitsArray *digarr=new AliTPCDigitsArray;
69digarr->Setup(param);
70cwd->cd();
71
72//loop over events
73Int_t nevmax=rl->GetNumberOfEvents();
74for(Int_t nev=0;nev<nevmax ;nev++){
75 rl->GetEvent(nev);
76 TTree* input_tree= tpcl->TreeD();//load tree with digits
77 digarr->ConnectTree(input_tree);
78 TTree *output_tree =tpcl->TreeR();//load output tree
79
80 AliTPCclustererKr *clusters = new AliTPCclustererKr();
81 clusters->SetParam(param);
82 clusters->SetInput(input_tree);
83 clusters->SetOutput(output_tree);
84 clusters->SetDigArr(digarr);
85
86//If you want to change the cluster finder parameters for MC there are
87//several of them:
88
89//1. signal threshold (everything below the given number is treated as 0)
90 clusters->SetMinAdc(3);
91
92//2. number of neighbouring timebins to be considered
93 clusters->SetMinTimeBins(2);
94
95//3. distance of the cluster center to the center of a pad in pad-padrow plane
96//(in cm). Remenber that this is still quantified by pad size.
97 clusters->SetMaxPadRangeCm(2.5);
98
99//4. distance of the cluster center to the center of a padrow in pad-padrow
100//plane (in cm). Remenber that this is still quantified by pad size.
101 clusters->SetMaxRowRangeCm(3.5);
102
103//5. distance of the cluster center to the max time bin on a pad (in tackts)
104//ie. fabs(centerT - time)<7
105 clusters->SetMaxTimeRange(7);
106
107//6. cut reduce peak at 0. There are noises which appear mostly as two
108//timebins on one pad.
109 clusters->SetValueToSize(3.5);
110
111
112 clusters->finderIO();
113 tpcl->WriteRecPoints("OVERWRITE");
114}
115delete rl;//cleans everything
116
117*
118********* DATA *********
119*
120
121To run clusterizaton for DATA for file named raw_data.root type:
122.x FindKrClustersRaw.C("raw_data.root")
123
124If you want to change some criteria do the following:
125
126//
127// remove Altro warnings
128//
129AliLog::SetClassDebugLevel("AliRawReaderDate",-5);
130AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5);
131AliLog::SetModuleDebugLevel("RAW",-5);
132
133//
134// Get database with noises
135//
136// char *ocdbpath = gSystem->Getenv("OCDB_PATH");
137char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB";
138if (ocdbpath==0){
139ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/";
140}
141AliCDBManager * man = AliCDBManager::Instance();
142man->SetDefaultStorage(ocdbpath);
143man->SetRun(0);
144AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
145AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
146
147//define tree
148TFile *hfile=new TFile("adc.root","RECREATE","ADC file");
149// Create a ROOT Tree
150TTree *mytree = new TTree("Kr","Krypton cluster tree");
151
152//define infput file
153const char *fileName="data.root";
154AliRawReader *reader = new AliRawReaderRoot(fileName);
155//AliRawReader *reader = new AliRawReaderDate(fileName);
156reader->Reset();
157AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader);
158stream->SelectRawData("TPC");
159
160//one general output
161AliTPCclustererKr *clusters = new AliTPCclustererKr();
162clusters->SetOutput(mytree);
163clusters->SetRecoParam(0);//standard reco parameters
164AliTPCParamSR *param=new AliTPCParamSR();
165clusters->SetParam(param);//TPC parameters(sectors, timebins, etc.)
166
167//set cluster finder parameters (from data):
168//1. zero suppression parameter
169 clusters->SetZeroSup(param->GetZeroSup());
170
171//2. first bin
172 clusters->SetFirstBin(60);
173
174//3. last bin
175 clusters->SetLastBin(950);
176
177//4. maximal noise
178 clusters->SetMaxNoiseAbs(2);
179
180//5. maximal amount of sigma of noise
181 clusters->SetMaxNoiseSigma(3);
182
183//The remaining parameters are the same paramters as for MC (see MC section
184//points 1-6)
185 clusters->SetMinAdc(3);
186 clusters->SetMinTimeBins(2);
187 clusters->SetMaxPadRangeCm(2.5);
188 clusters->SetMaxRowRangeCm(3.5);
189 clusters->SetMaxTimeRange(7);
190 clusters->SetValueToSize(3.5);
191
192while (reader->NextEvent()) {
193 clusters->FinderIO(reader);
194}
195
196hfile->Write();
197hfile->Close();
198delete stream;
199
200
201*/
202
203#include "AliTPCclustererKr.h"
204#include "AliTPCclusterKr.h"
205//#include <vector>
206#include <list>
207#include "TObject.h"
208#include "AliPadMax.h"
209#include "AliSimDigits.h"
210#include "AliTPCv4.h"
211#include "AliTPCParam.h"
212#include "AliTPCDigitsArray.h"
213#include "AliTPCvtpr.h"
214#include "AliTPCClustersRow.h"
215#include "TTree.h"
216#include "TH1F.h"
217#include "TH2F.h"
218#include "TTreeStream.h"
219
220#include "AliTPCTransform.h"
221
222//used in raw data finder
223#include "AliTPCROC.h"
224#include "AliTPCCalPad.h"
225#include "AliTPCAltroMapping.h"
226#include "AliTPCcalibDB.h"
227#include "AliTPCRawStreamV3.h"
228#include "AliTPCRecoParam.h"
229#include "AliTPCReconstructor.h"
230#include "AliRawReader.h"
231#include "AliTPCCalROC.h"
232#include "AliRawEventHeaderBase.h"
233
234using std::cerr;
235using std::cout;
236using std::endl;
237using std::list;
238ClassImp(AliTPCclustererKr)
239
240
241AliTPCclustererKr::AliTPCclustererKr()
242 :TObject(),
243 fRawData(kFALSE),
244 fInput(0),
245 fOutput(0),
246 fParam(0),
247 fDigarr(0),
248 fRecoParam(0),
249 fZeroSup(2),
250 fFirstBin(60),
251 fLastBin(950),
252 fMaxNoiseAbs(2),
253 fMaxNoiseSigma(3),
254 fMinAdc(3),
255 fMinTimeBins(2),
256// fMaxPadRange(4),
257// fMaxRowRange(3),
258 fMaxTimeRange(7),
259 fValueToSize(3.5),
260 fMaxPadRangeCm(2.5),
261 fMaxRowRangeCm(3.5),
262 fIsolCut(3),
263 fDebugLevel(-1),
264 fHistoRow(0),
265 fHistoPad(0),
266 fHistoTime(0),
267 fHistoRowPad(0),
268 fTimeStamp(0),
269 fRun(0)
270{
271//
272// default constructor
273//
274}
275
276AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr &param)
277 :TObject(),
278 fRawData(kFALSE),
279 fInput(0),
280 fOutput(0),
281 fParam(0),
282 fDigarr(0),
283 fRecoParam(0),
284 fZeroSup(2),
285 fFirstBin(60),
286 fLastBin(950),
287 fMaxNoiseAbs(2),
288 fMaxNoiseSigma(3),
289 fMinAdc(3),
290 fMinTimeBins(2),
291// fMaxPadRange(4),
292// fMaxRowRange(3),
293 fMaxTimeRange(7),
294 fValueToSize(3.5),
295 fMaxPadRangeCm(2.5),
296 fMaxRowRangeCm(3.5),
297 fIsolCut(3),
298 fDebugLevel(-1),
299 fHistoRow(0),
300 fHistoPad(0),
301 fHistoTime(0),
302 fHistoRowPad(0),
303 fTimeStamp(0),
304 fRun(0)
305{
306//
307// copy constructor
308//
309 fParam = param.fParam;
310 fRecoParam = param.fRecoParam;
311 fRawData = param.fRawData;
312 fInput = param.fInput ;
313 fOutput = param.fOutput;
314 fDigarr = param.fDigarr;
315 fZeroSup = param.fZeroSup ;
316 fFirstBin = param.fFirstBin ;
317 fLastBin = param.fLastBin ;
318 fMaxNoiseAbs = param.fMaxNoiseAbs ;
319 fMaxNoiseSigma = param.fMaxNoiseSigma ;
320 fMinAdc = param.fMinAdc;
321 fMinTimeBins = param.fMinTimeBins;
322// fMaxPadRange = param.fMaxPadRange ;
323// fMaxRowRange = param.fMaxRowRange ;
324 fMaxTimeRange = param.fMaxTimeRange;
325 fValueToSize = param.fValueToSize;
326 fMaxPadRangeCm = param.fMaxPadRangeCm;
327 fMaxRowRangeCm = param.fMaxRowRangeCm;
328 fIsolCut = param.fIsolCut;
329 fDebugLevel = param.fDebugLevel;
330 fHistoRow = param.fHistoRow ;
331 fHistoPad = param.fHistoPad ;
332 fHistoTime = param.fHistoTime;
333 fHistoRowPad = param.fHistoRowPad;
334 fTimeStamp = param.fTimeStamp;
335 fRun = param.fRun;
336
337}
338
339AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
340{
341 //
342 // assignment operator
343 //
344 if (this == &param) return (*this);
345
346 fParam = param.fParam;
347 fRecoParam = param.fRecoParam;
348 fRawData = param.fRawData;
349 fInput = param.fInput ;
350 fOutput = param.fOutput;
351 fDigarr = param.fDigarr;
352 fZeroSup = param.fZeroSup ;
353 fFirstBin = param.fFirstBin ;
354 fLastBin = param.fLastBin ;
355 fMaxNoiseAbs = param.fMaxNoiseAbs ;
356 fMaxNoiseSigma = param.fMaxNoiseSigma ;
357 fMinAdc = param.fMinAdc;
358 fMinTimeBins = param.fMinTimeBins;
359// fMaxPadRange = param.fMaxPadRange ;
360// fMaxRowRange = param.fMaxRowRange ;
361 fMaxTimeRange = param.fMaxTimeRange;
362 fValueToSize = param.fValueToSize;
363 fMaxPadRangeCm = param.fMaxPadRangeCm;
364 fMaxRowRangeCm = param.fMaxRowRangeCm;
365 fIsolCut = param.fIsolCut;
366 fDebugLevel = param.fDebugLevel;
367 fHistoRow = param.fHistoRow ;
368 fHistoPad = param.fHistoPad ;
369 fHistoTime = param.fHistoTime;
370 fHistoRowPad = param.fHistoRowPad;
371 fTimeStamp = param.fTimeStamp;
372 fRun = param.fRun;
373 return (*this);
374}
375
376AliTPCclustererKr::~AliTPCclustererKr()
377{
378 //
379 // destructor
380 //
381 delete fOutput;
382}
383
384void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
385{
386 //
387 // set reconstruction parameters
388 //
389 if (recoParam) {
390 fRecoParam = recoParam;
391 }else{
392 //set default parameters if not specified
393 fRecoParam = AliTPCReconstructor::GetRecoParam();
394 if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
395 }
396 return;
397}
398
399
400////____________________________________________________________________________
401//// I/O
402void AliTPCclustererKr::SetInput(TTree * tree)
403{
404 //
405 // set input tree with digits
406 //
407 fInput = tree;
408 if (!fInput->GetBranch("Segment")){
409 cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";
410 fInput=0;
411 return;
412 }
413}
414
415void AliTPCclustererKr::SetOutput(TTree * /*tree*/)
416{
417 //
418 //dummy
419 //
420 fOutput = new TTreeSRedirector("Krypton.root");
421}
422
423////____________________________________________________________________________
424//// with new I/O
425Int_t AliTPCclustererKr::FinderIO()
426{
427 // Krypton cluster finder for simulated events from MC
428
429 if (!fInput) {
430 Error("Digits2Clusters", "input tree not initialised");
431 return 10;
432 }
433
434 if (!fOutput) {
435 Error("Digits2Clusters", "output tree not initialised");
436 return 11;
437 }
438
439 FindClusterKrIO();
440 return 0;
441}
442
443
444
445Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)
446{
447 // Krypton cluster finder for the TPC raw data
448 // this method is unsing AliAltroRawStreamV3
449 // fParam must be defined before
450 if (!rawReader) return 1;
451 //
452 fRawData=kTRUE; //set flag to data
453
454 if (!fOutput) {
455 Error("Digits2Clusters", "output tree not initialised");
456 return 11;
457 }
458
459 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
460 // used later for memory allocation
461
462 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
463 if (eventHeader){
464 fTimeStamp = eventHeader->Get("Timestamp");
465 fRun = rawReader->GetRunNumber();
466 }
467
468
469 Bool_t isAltro=kFALSE;
470
471 AliTPCROC * roc = AliTPCROC::Instance();
472 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
473 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
474 //
475 AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
476
477 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
478 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
479 const Int_t kNS = kNIS + kNOS;//all sectors
480
481
482 //crate TPC view
483 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
484 digarr->Setup(fParam);//as usually parameters
485
486 for(Int_t iSec = 0; iSec < kNS; iSec++) {
487 AliTPCCalROC * noiseROC;
488 AliTPCCalROC noiseDummy(iSec);
489 if(noiseTPC==0x0){
490 noiseROC = &noiseDummy;//noise=0
491 }else{
492 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
493 }
494 Int_t nRows = 0; //number of rows in sector
495 Int_t nDDLs = 0; //number of DDLs
496 Int_t indexDDL = 0; //DDL index
497 if (iSec < kNIS) {
498 nRows = fParam->GetNRowLow();
499 nDDLs = 2;
500 indexDDL = iSec * 2;
501 }else {
502 nRows = fParam->GetNRowUp();
503 nDDLs = 4;
504 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
505 }
506
507 //
508 // Load the raw data for corresponding DDLs
509 //
510 rawReader->Reset();
511 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
512
513
514 while (input.NextDDL()){
515 // Allocate memory for rows in sector (pads(depends on row) x timebins)
516 if (!digarr->GetRow(iSec,0)){
517 for(Int_t iRow = 0; iRow < nRows; iRow++) {
518 digarr->CreateRow(iSec,iRow);
519 }//end loop over rows
520 }
521 //loop over pads
522 while ( input.NextChannel() ) {
523 Int_t iRow = input.GetRow();
524 Int_t iPad = input.GetPad();
525 //check row consistency
526 if (iRow < 0 ) continue;
527 if (iRow < 0 || iRow >= nRows){
528 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
529 iRow, 0, nRows -1));
530 continue;
531 }
532
533 //check pad consistency
534 if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
535 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
536 iPad, 0, roc->GetNPads(iSec,iRow) ));
537 continue;
538 }
539
540 //loop over bunches
541 while ( input.NextBunch() ){
542 Int_t startTbin = (Int_t)input.GetStartTimeBin();
543 Int_t bunchlength = (Int_t)input.GetBunchLength();
544 const UShort_t *sig = input.GetSignals();
545 isAltro=kTRUE;
546 for (Int_t iTime = 0; iTime<bunchlength; iTime++){
547 Int_t iTimeBin=startTbin-iTime;
548 //
549 if(fDebugLevel==72){
550 fHistoRow->Fill(iRow);
551 fHistoPad->Fill(iPad);
552 fHistoTime->Fill(iTimeBin);
553 fHistoRowPad->Fill(iPad,iRow);
554 }else if(fDebugLevel>=0&&fDebugLevel<72){
555 if(iSec==fDebugLevel){
556 fHistoRow->Fill(iRow);
557 fHistoPad->Fill(iPad);
558 fHistoTime->Fill(iTimeBin);
559 fHistoRowPad->Fill(iPad,iRow);
560 }
561 }else if(fDebugLevel==73){
562 if(iSec<36){
563 fHistoRow->Fill(iRow);
564 fHistoPad->Fill(iPad);
565 fHistoTime->Fill(iTimeBin);
566 fHistoRowPad->Fill(iPad,iRow);
567 }
568 }else if(fDebugLevel==74){
569 if(iSec>=36){
570 fHistoRow->Fill(iRow);
571 fHistoPad->Fill(iPad);
572 fHistoTime->Fill(iTimeBin);
573 fHistoRowPad->Fill(iPad,iRow);
574 }
575 }
576
577 //check time consistency
578 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
579 //cout<<iTimeBin<<endl;
580 continue;
581 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
582 iTimeBin, 0, fRecoParam->GetLastBin() -1));
583 }
584 //signal
585 Float_t signal=(Float_t)sig[iTime];
586 if (signal <= fZeroSup ||
587 iTimeBin < fFirstBin ||
588 iTimeBin > fLastBin
589 ) {
590 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
591 continue;
592 }
593 if (!noiseROC) continue;
594 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
595 if (noiseOnPad > fMaxNoiseAbs){
596 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
597 continue; // consider noisy pad as dead
598 }
599 if(signal <= fMaxNoiseSigma * noiseOnPad){
600 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
601 continue;
602 }
603 digarr->GetRow(iSec,iRow)->SetDigitFast(TMath::Nint(signal),iTimeBin,iPad);
604 }// end loop signals in bunch
605 }// end loop bunches
606 } // end loop pads
607 }// end ddl loop
608 }// end sector loop
609 SetDigArr(digarr);
610 if(isAltro) FindClusterKrIO();
611 delete digarr;
612
613 return 0;
614}
615
616void AliTPCclustererKr::CleanSector(Int_t sector){
617 //
618 // clean isolated digits
619 //
620 const Int_t kNRows=fParam->GetNRow(sector);//number of rows in sector
621 for(Int_t iRow=0; iRow<kNRows; ++iRow){
622 AliSimDigits *digrow;
623 if(fRawData){
624 digrow = (AliSimDigits*)fDigarr->GetRow(sector,iRow);//real data
625 }else{
626 digrow = (AliSimDigits*)fDigarr->LoadRow(sector,iRow);//MC
627 }
628 if(!digrow) continue;
629 digrow->ExpandBuffer(); //decrunch
630 const Int_t kNPads = digrow->GetNCols(); // number of pads
631 const Int_t kNTime = digrow->GetNRows(); // number of timebins
632 for(Int_t iPad=1;iPad<kNPads-1;iPad++){
633 Short_t* val = digrow->GetDigitsColumn(iPad);
634
635 for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
636 if (val[iTimeBin]<=0) continue;
637 if (val[iTimeBin-1]+val[iTimeBin+1]<fIsolCut) {val[iTimeBin]=0; continue;}
638 if (val[iTimeBin-kNTime]+val[iTimeBin+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
639 //
640 if (val[iTimeBin-1-kNTime]+val[iTimeBin+1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
641 if (val[iTimeBin+1-kNTime]+val[iTimeBin-1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
642
643 }
644 }
645 }
646}
647
648
649////____________________________________________________________________________
650Int_t AliTPCclustererKr::FindClusterKrIO()
651{
652
653 //
654 //fParam and fDigarr must be set to run this method
655 //
656
657 Int_t clusterCounter=0;
658 const Int_t nTotalSector=fParam->GetNSector();//number of sectors
659 for(Int_t iSec=0; iSec<nTotalSector; ++iSec){
660 CleanSector(iSec);
661
662 //vector of maxima for each sector
663 //std::vector<AliPadMax*> maximaInSector;
664 TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*
665
666 //
667 // looking for the maxima on the pad
668 //
669
670 const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector
671 for(Int_t iRow=0; iRow<kNRows; ++iRow){
672 AliSimDigits *digrow;
673 if(fRawData){
674 digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
675 }else{
676 digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
677 }
678 if(digrow){//if pointer exist
679 digrow->ExpandBuffer(); //decrunch
680 const Int_t kNPads = digrow->GetNCols(); // number of pads
681 const Int_t kNTime = digrow->GetNRows(); // number of timebins
682 for(Int_t iPad=0;iPad<kNPads;iPad++){
683
684 Int_t timeBinMax=-1;//timebin of maximum
685 Int_t valueMaximum=-1;//value of maximum in adc
686 Int_t increaseBegin=-1;//timebin when increase starts
687 Int_t sumAdc=0;//sum of adc on the pad in maximum surrounding
688 bool ifIncreaseBegin=true;//flag - check if increasing started
689 bool ifMaximum=false;//flag - check if it could be maximum
690 Short_t* val = digrow->GetDigitsColumn(iPad);
691 for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
692 if (!ifMaximum) {
693 if (val[iTimeBin]==-1) break; // 0 until the end
694 for( ; iTimeBin<kNTime-2&&val[iTimeBin]<fMinAdc ;iTimeBin++) {}
695 }
696 //
697 Short_t adc = val[iTimeBin];
698
699 if(adc<fMinAdc){//standard was 3 for fMinAdc
700 if(ifMaximum){
701 if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
702 timeBinMax=-1;
703 valueMaximum=-1;
704 increaseBegin=-1;
705 sumAdc=0;
706 ifIncreaseBegin=true;
707 ifMaximum=false;
708 continue;
709 }
710 //insert maximum, default values and set flags
711 //Double_t xCord,yCord;
712 //GetXY(iSec,iRow,iPad,xCord,yCord);
713 Double_t x[]={iRow,iPad,iTimeBin};
714 Int_t i[]={iSec};
715 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
716
717 transform->Transform(x,i,0,1);
718
719 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
720 timeBinMax,
721 iPad,
722 iRow,
723 x[0],//xCord,
724 x[1],//yCord,
725 x[2]/*timeBinMax*/),
726 increaseBegin,
727 iTimeBin-1,
728 sumAdc);
729 maximaInSector->AddLast(oneMaximum);
730
731 timeBinMax=-1;
732 valueMaximum=-1;
733 increaseBegin=-1;
734 sumAdc=0;
735 ifIncreaseBegin=true;
736 ifMaximum=false;
737 }
738 continue;
739 }
740
741
742
743
744
745
746 if(ifIncreaseBegin){
747 ifIncreaseBegin=false;
748 increaseBegin=iTimeBin;
749 }
750
751 if(adc>valueMaximum){
752 timeBinMax=iTimeBin;
753 valueMaximum=adc;
754 ifMaximum=true;
755 }
756 sumAdc+=adc;
757 if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
758 //at least 3 timebins
759 //insert maximum, default values and set flags
760 //Double_t xCord,yCord;
761 //GetXY(iSec,iRow,iPad,xCord,yCord);
762 Double_t x[]={iRow,iPad,iTimeBin};
763 Int_t i[]={iSec};
764 //AliTPCTransform trafo;
765 //trafo.Transform(x,i,0,1);
766
767 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
768
769 transform->Transform(x,i,0,1);
770
771 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
772 timeBinMax,
773 iPad,
774 iRow,
775 x[0],//xCord,
776 x[1],//yCord,
777 x[2]/*timeBinMax*/),
778 increaseBegin,
779 iTimeBin-1,
780 sumAdc);
781 maximaInSector->AddLast(oneMaximum);
782
783 timeBinMax=-1;
784 valueMaximum=-1;
785 increaseBegin=-1;
786 sumAdc=0;
787 ifIncreaseBegin=true;
788 ifMaximum=false;
789 continue;
790 }
791
792 }//end loop over timebins
793 }//end loop over pads
794// }else{
795// cout<<"Pointer does not exist!!"<<endl;
796 }//end if poiner exists
797 }//end loop over rows
798
799 MakeClusters(maximaInSector,iSec,clusterCounter);
800 //
801 maximaInSector->SetOwner(kTRUE);
802 maximaInSector->Delete();
803 delete maximaInSector;
804 }//end sector for
805 cout<<"Number of clusters in event: "<<clusterCounter<<endl;
806 return 0;
807}
808
809void AliTPCclustererKr::MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter){
810 //
811 // Make clusters
812 //
813
814 Int_t maxDig=0;
815 Int_t maxSumAdc=0;
816 Int_t maxTimeBin=0;
817 Int_t maxPad=0;
818 Int_t maxRow=0;
819 Double_t maxX=0;
820 Double_t maxY=0;
821 Double_t maxT=0;
822 Int_t entriesArr = maximaInSector->GetEntriesFast();
823 for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) {
824
825 AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1);
826 if (!mp1) continue;
827 AliTPCclusterKr clusterKr;
828
829 Int_t nUsedPads=1;
830 Int_t clusterValue=0;
831 clusterValue+=(mp1)->GetSum();
832 list<Int_t> nUsedRows;
833 nUsedRows.push_back((mp1)->GetRow());
834
835 maxDig =(mp1)->GetAdc() ;
836 maxSumAdc =(mp1)->GetSum() ;
837 maxTimeBin =(mp1)->GetTime();
838 maxPad =(mp1)->GetPad() ;
839 maxRow =(mp1)->GetRow() ;
840 maxX =(mp1)->GetX();
841 maxY =(mp1)->GetY();
842 maxT =(mp1)->GetT();
843
844 AliSimDigits *digrowTmp;
845 if(fRawData){
846 digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());
847 }else{
848 digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());
849 }
850
851 digrowTmp->ExpandBuffer(); //decrunch
852
853 for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){
854 Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad());
855 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT());
856 clusterKr.AddDigitToCluster(vtpr);
857 }
858 clusterKr.SetCenter();//set centr of the cluster
859
860 for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) {
861 AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2);
862 if (!mp2) continue;
863 if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX()) > fMaxPadRangeCm) continue;
864 if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue;
865 if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) > fMaxTimeRange) continue;
866
867 {
868 clusterValue+=(mp2)->GetSum();
869
870 nUsedPads++;
871 nUsedRows.push_back((mp2)->GetRow());
872
873 AliSimDigits *digrowTmp1;
874 if(fRawData){
875 digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());
876 }else{
877 digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());
878 }
879 digrowTmp1->ExpandBuffer(); //decrunch
880
881 for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){
882 Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad());
883 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT());
884 clusterKr.AddDigitToCluster(vtpr);
885 }
886
887 clusterKr.SetCenter();//set center of the cluster
888
889 //which one is bigger
890 if( (mp2)->GetAdc() > maxDig ){
891 maxDig =(mp2)->GetAdc() ;
892 maxSumAdc =(mp2)->GetSum() ;
893 maxTimeBin =(mp2)->GetTime();
894 maxPad =(mp2)->GetPad() ;
895 maxRow =(mp2)->GetRow() ;
896 maxX =(mp2)->GetX() ;
897 maxY =(mp2)->GetY() ;
898 maxT =(mp2)->GetT() ;
899 } else if ( (mp2)->GetAdc() == maxDig ){
900 if( (mp2)->GetSum() > maxSumAdc){
901 maxDig =(mp2)->GetAdc() ;
902 maxSumAdc =(mp2)->GetSum() ;
903 maxTimeBin =(mp2)->GetTime();
904 maxPad =(mp2)->GetPad() ;
905 maxRow =(mp2)->GetRow() ;
906 maxX =(mp2)->GetX() ;
907 maxY =(mp2)->GetY() ;
908 maxT =(mp2)->GetT() ;
909 }
910 }
911 delete maximaInSector->RemoveAt(it2);
912 }
913 }//inside loop
914 delete maximaInSector->RemoveAt(it1);
915 clusterKr.SetSize();
916 //through out clusters on the edge and noise
917 //if(clusterValue/clusterKr.fCluster.size()<fValueToSize)continue;
918 if(clusterValue/(clusterKr.GetSize())<fValueToSize)continue;
919
920 clusterKr.SetADCcluster(clusterValue);
921 clusterKr.SetNPads(nUsedPads);
922 clusterKr.SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxT));
923 clusterKr.SetSec(iSec);
924 clusterKr.SetSize();
925
926 nUsedRows.sort();
927 nUsedRows.unique();
928 clusterKr.SetNRows(nUsedRows.size());
929 clusterKr.SetCenter();
930
931 clusterKr.SetRMS();//Set pad,row,timebin RMS
932 clusterKr.Set1D();//Set size in pads and timebins
933
934 clusterKr.SetTimeStamp(fTimeStamp);
935 clusterKr.SetRun(fRun);
936
937 clusterCounter++;
938
939
940 //save each cluster into file
941 if (fOutput){
942 (*fOutput)<<"Kr"<<
943 "Cl.="<<&clusterKr<<
944 "\n";
945 }
946 //end of save each cluster into file adc.root
947 }//outer loop
948}
949
950
951
952////____________________________________________________________________________
953
954
955void AliTPCclustererKr::GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob){
956 //
957 //gives global XY coordinate of the pad
958 //
959
960 Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
961
962 Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
963 Float_t padXSize;
964 if(sec<fParam->GetNInnerSector())padXSize=0.4;
965 else padXSize=0.6;
966 Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad
967
968 Float_t sin,cos;
969 fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
970
971 Double_t xGlob1 = xLocal * cos + yLocal * sin;
972 Double_t yGlob1 = -xLocal * sin + yLocal * cos;
973
974
975 Double_t rot=0;
976 rot=TMath::Pi()/2.;
977
978 xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);
979 yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);
980
981 yGlob=-1*yGlob;
982 if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;
983
984
985 return;
986}