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