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