Just something left from v2...now clean
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererKr.cxx
CommitLineData
87ef42d4 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
557ab431 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
87ef42d4 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("AliTPCRawStream",-5);
130AliLog::SetClassDebugLevel("AliRawReaderDate",-5);
131AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5);
132AliLog::SetModuleDebugLevel("RAW",-5);
133
134//
135// Get database with noises
136//
137// char *ocdbpath = gSystem->Getenv("OCDB_PATH");
138char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB";
139if (ocdbpath==0){
140ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/";
141}
142AliCDBManager * man = AliCDBManager::Instance();
143man->SetDefaultStorage(ocdbpath);
144man->SetRun(0);
145AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
146AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
147
148//define tree
149TFile *hfile=new TFile("adc.root","RECREATE","ADC file");
150// Create a ROOT Tree
151TTree *mytree = new TTree("Kr","Krypton cluster tree");
152
153//define infput file
154const char *fileName="data.root";
155AliRawReader *reader = new AliRawReaderRoot(fileName);
156//AliRawReader *reader = new AliRawReaderDate(fileName);
157reader->Reset();
158AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader);
159stream->SelectRawData("TPC");
160
161//one general output
162AliTPCclustererKr *clusters = new AliTPCclustererKr();
163clusters->SetOutput(mytree);
164clusters->SetRecoParam(0);//standard reco parameters
165AliTPCParamSR *param=new AliTPCParamSR();
166clusters->SetParam(param);//TPC parameters(sectors, timebins, etc.)
167
168//set cluster finder parameters (from data):
169//1. zero suppression parameter
170 clusters->SetZeroSup(param->GetZeroSup());
171
172//2. first bin
173 clusters->SetFirstBin(60);
174
175//3. last bin
176 clusters->SetLastBin(950);
177
178//4. maximal noise
179 clusters->SetMaxNoiseAbs(2);
180
181//5. maximal amount of sigma of noise
182 clusters->SetMaxNoiseSigma(3);
183
184//The remaining parameters are the same paramters as for MC (see MC section
185//points 1-6)
186 clusters->SetMinAdc(3);
187 clusters->SetMinTimeBins(2);
188 clusters->SetMaxPadRangeCm(2.5);
189 clusters->SetMaxRowRangeCm(3.5);
190 clusters->SetMaxTimeRange(7);
191 clusters->SetValueToSize(3.5);
192
193while (reader->NextEvent()) {
194 clusters->FinderIO(reader);
195}
196
197hfile->Write();
198hfile->Close();
199delete stream;
200
201
202*/
203
204#include "AliTPCclustererKr.h"
205#include "AliTPCclusterKr.h"
206//#include <vector>
207#include <list>
208#include "TObject.h"
209#include "AliPadMax.h"
210#include "AliSimDigits.h"
211#include "AliTPCv4.h"
212#include "AliTPCParam.h"
213#include "AliTPCDigitsArray.h"
214#include "AliTPCvtpr.h"
215#include "AliTPCClustersRow.h"
216#include "TTree.h"
217#include "TH1F.h"
218#include "TH2F.h"
26a97212 219#include "TTreeStream.h"
87ef42d4 220
dd53103b 221#include "AliTPCTransform.h"
222
87ef42d4 223//used in raw data finder
224#include "AliTPCROC.h"
225#include "AliTPCCalPad.h"
226#include "AliTPCAltroMapping.h"
227#include "AliTPCcalibDB.h"
228#include "AliTPCRawStream.h"
229#include "AliTPCRecoParam.h"
230#include "AliTPCReconstructor.h"
231#include "AliRawReader.h"
232#include "AliTPCCalROC.h"
233
234ClassImp(AliTPCclustererKr)
235
236
237AliTPCclustererKr::AliTPCclustererKr()
238 :TObject(),
239 fRawData(kFALSE),
87ef42d4 240 fInput(0),
241 fOutput(0),
242 fParam(0),
243 fDigarr(0),
244 fRecoParam(0),
245 fZeroSup(2),
246 fFirstBin(60),
247 fLastBin(950),
248 fMaxNoiseAbs(2),
249 fMaxNoiseSigma(3),
250 fMinAdc(3),
251 fMinTimeBins(2),
252// fMaxPadRange(4),
253// fMaxRowRange(3),
254 fMaxTimeRange(7),
255 fValueToSize(3.5),
256 fMaxPadRangeCm(2.5),
257 fMaxRowRangeCm(3.5),
258 fDebugLevel(-1),
259 fHistoRow(0),
260 fHistoPad(0),
261 fHistoTime(0),
262 fHistoRowPad(0)
263{
264//
265// default constructor
266//
267}
268
269AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr &param)
270 :TObject(),
271 fRawData(kFALSE),
87ef42d4 272 fInput(0),
273 fOutput(0),
274 fParam(0),
275 fDigarr(0),
276 fRecoParam(0),
277 fZeroSup(2),
278 fFirstBin(60),
279 fLastBin(950),
280 fMaxNoiseAbs(2),
281 fMaxNoiseSigma(3),
282 fMinAdc(3),
283 fMinTimeBins(2),
284// fMaxPadRange(4),
285// fMaxRowRange(3),
286 fMaxTimeRange(7),
287 fValueToSize(3.5),
288 fMaxPadRangeCm(2.5),
289 fMaxRowRangeCm(3.5),
290 fDebugLevel(-1),
291 fHistoRow(0),
292 fHistoPad(0),
293 fHistoTime(0),
294 fHistoRowPad(0)
295{
296//
297// copy constructor
298//
299 fParam = param.fParam;
300 fRecoParam = param.fRecoParam;
301 fRawData = param.fRawData;
87ef42d4 302 fInput = param.fInput ;
303 fOutput = param.fOutput;
304 fDigarr = param.fDigarr;
305 fZeroSup = param.fZeroSup ;
306 fFirstBin = param.fFirstBin ;
307 fLastBin = param.fLastBin ;
308 fMaxNoiseAbs = param.fMaxNoiseAbs ;
309 fMaxNoiseSigma = param.fMaxNoiseSigma ;
310 fMinAdc = param.fMinAdc;
311 fMinTimeBins = param.fMinTimeBins;
312// fMaxPadRange = param.fMaxPadRange ;
313// fMaxRowRange = param.fMaxRowRange ;
314 fMaxTimeRange = param.fMaxTimeRange;
315 fValueToSize = param.fValueToSize;
316 fMaxPadRangeCm = param.fMaxPadRangeCm;
317 fMaxRowRangeCm = param.fMaxRowRangeCm;
318 fDebugLevel = param.fDebugLevel;
319 fHistoRow = param.fHistoRow ;
320 fHistoPad = param.fHistoPad ;
321 fHistoTime = param.fHistoTime;
322 fHistoRowPad = param.fHistoRowPad;
323
324}
325
326AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
327{
328 //
329 // assignment operator
330 //
331 fParam = param.fParam;
332 fRecoParam = param.fRecoParam;
333 fRawData = param.fRawData;
87ef42d4 334 fInput = param.fInput ;
335 fOutput = param.fOutput;
336 fDigarr = param.fDigarr;
337 fZeroSup = param.fZeroSup ;
338 fFirstBin = param.fFirstBin ;
339 fLastBin = param.fLastBin ;
340 fMaxNoiseAbs = param.fMaxNoiseAbs ;
341 fMaxNoiseSigma = param.fMaxNoiseSigma ;
342 fMinAdc = param.fMinAdc;
343 fMinTimeBins = param.fMinTimeBins;
344// fMaxPadRange = param.fMaxPadRange ;
345// fMaxRowRange = param.fMaxRowRange ;
346 fMaxTimeRange = param.fMaxTimeRange;
347 fValueToSize = param.fValueToSize;
348 fMaxPadRangeCm = param.fMaxPadRangeCm;
349 fMaxRowRangeCm = param.fMaxRowRangeCm;
350 fDebugLevel = param.fDebugLevel;
351 fHistoRow = param.fHistoRow ;
352 fHistoPad = param.fHistoPad ;
353 fHistoTime = param.fHistoTime;
354 fHistoRowPad = param.fHistoRowPad;
355 return (*this);
356}
357
358AliTPCclustererKr::~AliTPCclustererKr()
359{
360 //
361 // destructor
362 //
26a97212 363 delete fOutput;
87ef42d4 364}
365
366void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
367{
368 //
369 // set reconstruction parameters
370 //
371 if (recoParam) {
372 fRecoParam = recoParam;
373 }else{
374 //set default parameters if not specified
375 fRecoParam = AliTPCReconstructor::GetRecoParam();
376 if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
377 }
378 return;
379}
380
381
382////____________________________________________________________________________
383//// I/O
384void AliTPCclustererKr::SetInput(TTree * tree)
385{
386 //
387 // set input tree with digits
388 //
389 fInput = tree;
390 if (!fInput->GetBranch("Segment")){
391 cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";
392 fInput=0;
393 return;
394 }
395}
396
26a97212 397void AliTPCclustererKr::SetOutput(TTree * /*tree*/)
87ef42d4 398{
399 //
26a97212 400 //dummy
87ef42d4 401 //
26a97212 402 fOutput = new TTreeSRedirector("Krypton.root");
87ef42d4 403}
404
405////____________________________________________________________________________
406//// with new I/O
407Int_t AliTPCclustererKr::FinderIO()
408{
409 // Krypton cluster finder for simulated events from MC
410
411 if (!fInput) {
412 Error("Digits2Clusters", "input tree not initialised");
413 return 10;
414 }
415
416 if (!fOutput) {
417 Error("Digits2Clusters", "output tree not initialised");
418 return 11;
419 }
420
421 FindClusterKrIO();
422 return 0;
423}
424
425Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)
426{
427 // Krypton cluster finder for the TPC raw data
428 //
429 // fParam must be defined before
430
431 if(rawReader)fRawData=kTRUE; //set flag to data
432
433 if (!fOutput) {
434 Error("Digits2Clusters", "output tree not initialised");
435 return 11;
436 }
437
438 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
439 // used later for memory allocation
440
441 Bool_t isAltro=kFALSE;
442
443 AliTPCROC * roc = AliTPCROC::Instance();
444 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
445 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
446 //
447 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
448
449 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
450 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
451 const Int_t kNS = kNIS + kNOS;//all sectors
452
87ef42d4 453
454 //crate TPC view
455 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
456 digarr->Setup(fParam);//as usually parameters
457
458 //
459 // Loop over sectors
460 //
461 for(Int_t iSec = 0; iSec < kNS; iSec++) {
462 AliTPCCalROC * noiseROC;
463 if(noiseTPC==0x0){
464 noiseROC =new AliTPCCalROC(iSec);//noise=0
465 }
466 else{
467 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
468 }
469 Int_t nRows = 0; //number of rows in sector
470 Int_t nDDLs = 0; //number of DDLs
471 Int_t indexDDL = 0; //DDL index
472 if (iSec < kNIS) {
473 nRows = fParam->GetNRowLow();
474 nDDLs = 2;
475 indexDDL = iSec * 2;
476 }else {
477 nRows = fParam->GetNRowUp();
478 nDDLs = 4;
479 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
480 }
481
482 //
483 // Load the raw data for corresponding DDLs
484 //
485 rawReader->Reset();
486 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
487
488 if(input.Next()) {
489 isAltro=kTRUE;
490 // Allocate memory for rows in sector (pads(depends on row) x timebins)
491 for(Int_t iRow = 0; iRow < nRows; iRow++) {
557ab431 492 digarr->CreateRow(iSec,iRow);
87ef42d4 493 }//end loop over rows
494 }
495 rawReader->Reset();
496 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
497
498 //
499 // Begin loop over altro data
500 //
501 while (input.Next()) {
557ab431 502
87ef42d4 503 //check sector consistency
504 if (input.GetSector() != iSec)
505 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSec,input.GetSector()));
506
557ab431 507 Int_t iRow = input.GetRow();
508 Int_t iPad = input.GetPad();
509 Int_t iTimeBin = input.GetTime();
87ef42d4 510
557ab431 511 //
87ef42d4 512 if(fDebugLevel==72){
513 fHistoRow->Fill(iRow);
514 fHistoPad->Fill(iPad);
515 fHistoTime->Fill(iTimeBin);
516 fHistoRowPad->Fill(iPad,iRow);
517 }else if(fDebugLevel>=0&&fDebugLevel<72){
518 if(iSec==fDebugLevel){
519 fHistoRow->Fill(iRow);
520 fHistoPad->Fill(iPad);
521 fHistoTime->Fill(iTimeBin);
522 fHistoRowPad->Fill(iPad,iRow);
523 }
524 }else if(fDebugLevel==73){
525 if(iSec<36){
526 fHistoRow->Fill(iRow);
527 fHistoPad->Fill(iPad);
528 fHistoTime->Fill(iTimeBin);
529 fHistoRowPad->Fill(iPad,iRow);
530 }
531 }else if(fDebugLevel==74){
532 if(iSec>=36){
533 fHistoRow->Fill(iRow);
534 fHistoPad->Fill(iPad);
535 fHistoTime->Fill(iTimeBin);
536 fHistoRowPad->Fill(iPad,iRow);
537 }
538 }
539
540 //check row consistency
26a97212 541 if (iRow < 0 ) continue;
87ef42d4 542 if (iRow < 0 || iRow >= nRows){
543 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
544 iRow, 0, nRows -1));
545 continue;
546 }
547
548 //check pad consistency
557ab431 549 if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
87ef42d4 550 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
551 iPad, 0, roc->GetNPads(iSec,iRow) ));
552 continue;
553 }
554
555 //check time consistency
556 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
557 //cout<<iTimeBin<<endl;
558 continue;
559 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
560 iTimeBin, 0, fRecoParam->GetLastBin() -1));
561 }
562
563 //signal
564 Int_t signal = input.GetSignal();
565 if (signal <= fZeroSup ||
566 iTimeBin < fFirstBin ||
567 iTimeBin > fLastBin
568 ) {
569 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
570 continue;
571 }
572 if (!noiseROC) continue;
573 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
574 if (noiseOnPad > fMaxNoiseAbs){
575 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
576 continue; // consider noisy pad as dead
577 }
578 if(signal <= fMaxNoiseSigma * noiseOnPad){
579 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
580 continue;
581 }
557ab431 582 digarr->GetRow(iSec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
87ef42d4 583 }//end of loop over altro data
584 }//end of loop over sectors
585
586 SetDigArr(digarr);
587 if(isAltro) FindClusterKrIO();
588 delete digarr;
589
590 return 0;
591}
592
593////____________________________________________________________________________
594Int_t AliTPCclustererKr::FindClusterKrIO()
595{
87ef42d4 596
557ab431 597 //
598 //fParam and fDigarr must be set to run this method
599 //
87ef42d4 600
601 Int_t clusterCounter=0;
557ab431 602 const Int_t nTotalSector=fParam->GetNSector();//number of sectors
603 for(Int_t iSec=0; iSec<nTotalSector; ++iSec){
87ef42d4 604
605 //vector of maxima for each sector
606 //std::vector<AliPadMax*> maximaInSector;
607 TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*
608
609 //
610 // looking for the maxima on the pad
611 //
612
557ab431 613 const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector
614 for(Int_t iRow=0; iRow<kNRows; ++iRow){
87ef42d4 615 AliSimDigits *digrow;
616 if(fRawData){
617 digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
618 }else{
619 digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
620 }
621 if(digrow){//if pointer exist
622 digrow->ExpandBuffer(); //decrunch
557ab431 623 const Int_t kNPads = digrow->GetNCols(); // number of pads
624 const Int_t kNTime = digrow->GetNRows(); // number of timebins
625 for(Int_t iPad=0;iPad<kNPads;iPad++){
87ef42d4 626
557ab431 627 Int_t timeBinMax=-1;//timebin of maximum
628 Int_t valueMaximum=-1;//value of maximum in adc
629 Int_t increaseBegin=-1;//timebin when increase starts
630 Int_t sumAdc=0;//sum of adc on the pad in maximum surrounding
87ef42d4 631 bool ifIncreaseBegin=true;//flag - check if increasing started
632 bool ifMaximum=false;//flag - check if it could be maximum
557ab431 633 Short_t* val = digrow->GetDigitsColumn(iPad);
634 for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
635 if (!ifMaximum) {
636 if (val[iTimeBin]==-1) break; // 0 until the end
637 for( ; iTimeBin<kNTime-2&&val[iTimeBin]<fMinAdc ;iTimeBin++);
638 }
639 //
640 Short_t adc = val[iTimeBin];
87ef42d4 641 if(adc<fMinAdc){//standard was 3
642 if(ifMaximum){
643 if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
644 timeBinMax=-1;
645 valueMaximum=-1;
646 increaseBegin=-1;
647 sumAdc=0;
648 ifIncreaseBegin=true;
649 ifMaximum=false;
650 continue;
651 }
652 //insert maximum, default values and set flags
dd53103b 653 //Double_t xCord,yCord;
654 //GetXY(iSec,iRow,iPad,xCord,yCord);
655 Double_t x[]={iRow,iPad,iTimeBin};
656 Int_t i[]={iSec};
657 AliTPCTransform trafo;
658 trafo.Transform(x,i,0,1);
557ab431 659
87ef42d4 660 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
661 timeBinMax,
662 iPad,
663 iRow,
dd53103b 664 x[0],//xCord,
665 x[1],//yCord,
666 x[2]/*timeBinMax*/),
87ef42d4 667 increaseBegin,
668 iTimeBin-1,
669 sumAdc);
87ef42d4 670 maximaInSector->AddLast(oneMaximum);
671
672 timeBinMax=-1;
673 valueMaximum=-1;
674 increaseBegin=-1;
675 sumAdc=0;
676 ifIncreaseBegin=true;
677 ifMaximum=false;
678 }
679 continue;
680 }
681
682 if(ifIncreaseBegin){
683 ifIncreaseBegin=false;
684 increaseBegin=iTimeBin;
685 }
686
687 if(adc>valueMaximum){
688 timeBinMax=iTimeBin;
689 valueMaximum=adc;
690 ifMaximum=true;
691 }
692 sumAdc+=adc;
693 if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
694 //at least 3 timebins
695 //insert maximum, default values and set flags
dd53103b 696 //Double_t xCord,yCord;
697 //GetXY(iSec,iRow,iPad,xCord,yCord);
698 Double_t x[]={iRow,iPad,iTimeBin};
699 Int_t i[]={iSec};
700 AliTPCTransform trafo;
701 trafo.Transform(x,i,0,1);
87ef42d4 702 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
703 timeBinMax,
704 iPad,
705 iRow,
dd53103b 706 x[0],//xCord,
707 x[1],//yCord,
708 x[2]/*timeBinMax*/),
87ef42d4 709 increaseBegin,
710 iTimeBin-1,
711 sumAdc);
87ef42d4 712 maximaInSector->AddLast(oneMaximum);
713
714 timeBinMax=-1;
715 valueMaximum=-1;
716 increaseBegin=-1;
717 sumAdc=0;
718 ifIncreaseBegin=true;
719 ifMaximum=false;
720 continue;
721 }
722
723 }//end loop over timebins
724 }//end loop over pads
725// }else{
726// cout<<"Pointer does not exist!!"<<endl;
727 }//end if poiner exists
728 }//end loop over rows
729
557ab431 730 MakeClusters(maximaInSector,iSec,clusterCounter);
87ef42d4 731 //
557ab431 732 maximaInSector->SetOwner(kTRUE);
733 maximaInSector->Delete();
734 delete maximaInSector;
735 }//end sector for
736 cout<<"Number of clusters in event: "<<clusterCounter<<endl;
737 return 0;
738}
87ef42d4 739
557ab431 740void AliTPCclustererKr::MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter){
741 //
742 // Make clusters
743 //
87ef42d4 744
557ab431 745 Int_t maxDig=0;
746 Int_t maxSumAdc=0;
747 Int_t maxTimeBin=0;
748 Int_t maxPad=0;
749 Int_t maxRow=0;
750 Double_t maxX=0;
751 Double_t maxY=0;
752 Double_t maxT=0;
753 Int_t entriesArr = maximaInSector->GetEntriesFast();
754 for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) {
755
756 AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1);
757 if (!mp1) continue;
758 AliTPCclusterKr clusterKr;
759
760 Int_t nUsedPads=1;
761 Int_t clusterValue=0;
762 clusterValue+=(mp1)->GetSum();
763 list<Int_t> nUsedRows;
764 nUsedRows.push_back((mp1)->GetRow());
765
766 maxDig =(mp1)->GetAdc() ;
767 maxSumAdc =(mp1)->GetSum() ;
768 maxTimeBin =(mp1)->GetTime();
769 maxPad =(mp1)->GetPad() ;
770 maxRow =(mp1)->GetRow() ;
771 maxX =(mp1)->GetX();
772 maxY =(mp1)->GetY();
773 maxT =(mp1)->GetT();
774
775 AliSimDigits *digrowTmp;
776 if(fRawData){
777 digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());
778 }else{
779 digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());
780 }
781
782 digrowTmp->ExpandBuffer(); //decrunch
783
784 for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){
785 Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad());
786 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT());
787 clusterKr.AddDigitToCluster(vtpr);
788 }
789 clusterKr.SetCenter();//set centr of the cluster
790
791 for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) {
792 AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2);
793 if (!mp2) continue;
794 if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX())>fMaxPadRangeCm) continue;
795 if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue;
796 if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) >fMaxTimeRange) continue;
797
798 {
799 clusterValue+=(mp2)->GetSum();
800
801 nUsedPads++;
802 nUsedRows.push_back((mp2)->GetRow());
803
804 AliSimDigits *digrowTmp1;
805 if(fRawData){
806 digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());
807 }else{
808 digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());
809 }
810 digrowTmp1->ExpandBuffer(); //decrunch
811
812 for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){
813 Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad());
814 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT());
815 clusterKr.AddDigitToCluster(vtpr);
816 }
817
818 clusterKr.SetCenter();//set center of the cluster
819
820 //which one is bigger
821 if( (mp2)->GetAdc() > maxDig ){
822 maxDig =(mp2)->GetAdc() ;
823 maxSumAdc =(mp2)->GetSum() ;
824 maxTimeBin =(mp2)->GetTime();
825 maxPad =(mp2)->GetPad() ;
826 maxRow =(mp2)->GetRow() ;
827 maxX =(mp2)->GetX() ;
828 maxY =(mp2)->GetY() ;
829 maxT =(mp2)->GetT() ;
830 } else if ( (mp2)->GetAdc() == maxDig ){
831 if( (mp2)->GetSum() > maxSumAdc){
87ef42d4 832 maxDig =(mp2)->GetAdc() ;
833 maxSumAdc =(mp2)->GetSum() ;
834 maxTimeBin =(mp2)->GetTime();
835 maxPad =(mp2)->GetPad() ;
836 maxRow =(mp2)->GetRow() ;
837 maxX =(mp2)->GetX() ;
838 maxY =(mp2)->GetY() ;
dd53103b 839 maxT =(mp2)->GetT() ;
87ef42d4 840 }
87ef42d4 841 }
557ab431 842 delete maximaInSector->RemoveAt(it2);
26a97212 843 }
557ab431 844 }//inside loop
845 delete maximaInSector->RemoveAt(it1);
846 clusterKr.SetSize();
847 //through out clusters on the edge and noise
848 //if(clusterValue/clusterKr.fCluster.size()<fValueToSize)continue;
849 if(clusterValue/(clusterKr.GetSize())<fValueToSize)continue;
850
851 clusterKr.SetADCcluster(clusterValue);
852 clusterKr.SetNPads(nUsedPads);
853 clusterKr.SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxT));
854 clusterKr.SetSec(iSec);
855 clusterKr.SetSize();
856
857 nUsedRows.sort();
858 nUsedRows.unique();
859 clusterKr.SetNRows(nUsedRows.size());
860 clusterKr.SetCenter();
861
862 clusterCounter++;
863
864
865 //save each cluster into file
866 if (fOutput){
867 (*fOutput)<<"Kr"<<
868 "Cl.="<<&clusterKr<<
869 "\n";
870 }
871 //end of save each cluster into file adc.root
872 }//outer loop
87ef42d4 873}
874
557ab431 875
876
87ef42d4 877////____________________________________________________________________________
878
879
557ab431 880void AliTPCclustererKr::GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob){
87ef42d4 881 //
882 //gives global XY coordinate of the pad
883 //
884
885 Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
886
557ab431 887 Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
87ef42d4 888 Float_t padXSize;
889 if(sec<fParam->GetNInnerSector())padXSize=0.4;
890 else padXSize=0.6;
891 Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad
892
893 Float_t sin,cos;
894 fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
895
896 Double_t xGlob1 = xLocal * cos + yLocal * sin;
897 Double_t yGlob1 = -xLocal * sin + yLocal * cos;
898
899
900 Double_t rot=0;
901 rot=TMath::Pi()/2.;
902
903 xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);
904 yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);
905
906 yGlob=-1*yGlob;
907 if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;
908
909
910 return;
911}