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