TPC/AliTPCclustererKr.cxx

   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
  25 Instruction - how to use that:\r
  26 There are two macros prepared. One is for preparing clusters from MC \r
  27 samples:  FindKrClusters.C. The output is kept in TPC.RecPoints.root.\r
  28 The other macro is prepared for data analysis: FindKrClustersRaw.C. \r
  29 The output is created for each processed file in root file named adc.root. \r
  30 For each data subsample the same named file is created. So be careful \r
  31 do not overwrite them. \r
  32 \r
  33 *\r
  34 **** MC ****\r
  35 *\r
  36 \r
  37 To run clusterizaton for MC type:\r
  38 .x FindKrClusters.C\r
  39 \r
  40 If you don't want to use the standard selection criteria then you \r
  41 have to do following:\r
  42 \r
  43 // load the standard setup\r
  44 AliRunLoader* rl = AliRunLoader::Open("galice.root");\r
  45 AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");\r
  46 tpcl->LoadDigits();\r
  47 rl->LoadgAlice();\r
  48 gAlice=rl->GetAliRun();\r
  49 TDirectory *cwd = gDirectory;\r
  50 AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC");\r
  51 Int_t ver = tpc->IsVersion();\r
  52 rl->CdGAFile();\r
  53 AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");\r
  54 AliTPCDigitsArray *digarr=new AliTPCDigitsArray;\r
  55 digarr->Setup(param);\r
  56 cwd->cd();\r
  57 \r
  58 //loop over events\r
  59 Int_t nevmax=rl->GetNumberOfEvents();\r
  60 for(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
 101 delete rl;//cleans everything\r
 102 \r
 103 *\r
 104 ********* DATA *********\r
 105 *\r
 106 \r
 107 To run clusterizaton for DATA for file named raw_data.root type:\r
 108 .x FindKrClustersRaw.C("raw_data.root")\r
 109 \r
 110 If you want to change some criteria do the following:\r
 111 \r
 112 //\r
 113 // remove Altro warnings\r
 114 //\r
 115 AliLog::SetClassDebugLevel("AliTPCRawStream",-5);\r
 116 AliLog::SetClassDebugLevel("AliRawReaderDate",-5);\r
 117 AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5);\r
 118 AliLog::SetModuleDebugLevel("RAW",-5);\r
 119 \r
 120 //\r
 121 // Get database with noises\r
 122 //\r
 123 //  char *ocdbpath = gSystem->Getenv("OCDB_PATH");\r
 124 char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB";\r
 125 if (ocdbpath==0){\r
 126 ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/";\r
 127 }\r
 128 AliCDBManager * man = AliCDBManager::Instance();\r
 129 man->SetDefaultStorage(ocdbpath);\r
 130 man->SetRun(0);\r
 131 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();\r
 132 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();\r
 133 \r
 134 //define tree\r
 135 TFile *hfile=new TFile("adc.root","RECREATE","ADC file");\r
 136 // Create a ROOT Tree\r
 137 TTree *mytree = new TTree("Kr","Krypton cluster tree");\r
 138 \r
 139 //define infput file\r
 140 const char *fileName="data.root";\r
 141 AliRawReader *reader = new AliRawReaderRoot(fileName);\r
 142 //AliRawReader *reader = new AliRawReaderDate(fileName);\r
 143 reader->Reset();\r
 144 AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader);\r
 145 stream->SelectRawData("TPC");\r
 146 \r
 147 //one general output\r
 148 AliTPCclustererKr *clusters = new AliTPCclustererKr();\r
 149 clusters->SetOutput(mytree);\r
 150 clusters->SetRecoParam(0);//standard reco parameters\r
 151 AliTPCParamSR *param=new AliTPCParamSR();\r
 152 clusters->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
 179 while (reader->NextEvent()) {\r
 180   clusters->FinderIO(reader);\r
 181 }\r
 182 \r
 183 hfile->Write();\r
 184 hfile->Close();\r
 185 delete stream;\r
 186 \r
 187 \r
 188 */\r
 189 \r
 190 #include "AliTPCclustererKr.h"\r
 191 #include "AliTPCclusterKr.h"\r
 192 //#include <vector>\r
 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
 203 #include "TH1F.h"\r
 204 #include "TH2F.h"\r
 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
 217 ClassImp(AliTPCclustererKr)\r
 218 \r
 219 \r
 220 AliTPCclustererKr::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
 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
 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
 247 {\r
 248 //\r
 249 // default constructor\r
 250 //\r
 251 }\r
 252 \r
 253 AliTPCclustererKr::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
 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
 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
 280 {\r
 281 //\r
 282 // copy constructor\r
 283 //\r
 284   fParam = param.fParam;\r
 285   fRecoParam = param.fRecoParam;\r
 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
 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
 310 } \r
 311 \r
 312 AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)\r
 313 {\r
 314   //\r
 315   // assignment operator\r
 316   //\r
 317   fParam = param.fParam;\r
 318   fRecoParam = param.fRecoParam;\r
 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
 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
 342   return (*this);\r
 343 }\r
 344 \r
 345 AliTPCclustererKr::~AliTPCclustererKr()\r
 346 {\r
 347   //\r
 348   // destructor\r
 349   //\r
 350 }\r
 351 \r
 352 void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)\r
 353 {\r
 354   //\r
 355   // set reconstruction parameters\r
 356   //\r
 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
 370 void 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
 383 void 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
 398 Int_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
 416 Int_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
 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
 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
 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
 504       Short_t iRow = input.GetRow();\r
 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
 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
 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
 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
 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
 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
 591 Int_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
 610     //std::vector<AliPadMax*> maximaInSector;\r
 611     TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*\r
 612 \r
 613     //\r
 614     //  looking for the maxima on the pad\r
 615     //\r
 616 \r
 617     const Short_t kNRows=fParam->GetNRow(iSec);//number of rows in sector\r
 618     for(Short_t iRow=0; iRow<kNRows; ++iRow){\r
 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
 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
 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
 638           for(Short_t iTimeBin=0;iTimeBin<kNTime;iTimeBin++){\r
 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
 664                 //maximaInSector.push_back(oneMaximum);\r
 665                 maximaInSector->AddLast(oneMaximum);\r
 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
 688             if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge\r
 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
 703               //maximaInSector.push_back(oneMaximum);\r
 704               maximaInSector->AddLast(oneMaximum);\r
 705                 \r
 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
 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
 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
 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
 744       AliTPCclusterKr *tmp=new AliTPCclusterKr();\r
 745       \r
 746       Short_t nUsedPads=1;\r
 747       Int_t clusterValue=0;\r
 748       clusterValue+=(mp1)->GetSum();\r
 749       list<Short_t> nUsedRows;\r
 750       nUsedRows.push_back((mp1)->GetRow());\r
 751 \r
 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
 759 \r
 760       AliSimDigits *digrowTmp;\r
 761       if(fRawData){\r
 762         digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());\r
 763       }else{\r
 764         digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());\r
 765       }\r
 766 \r
 767       digrowTmp->ExpandBuffer(); //decrunch\r
 768 \r
 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
 774       }\r
 775       tmp->SetCenter();//set centr of the cluster\r
 776       \r
 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
 790 //      if(abs(maxRow - (*mp2)->GetRow()) < fMaxRowRange && //3\r
 791 //         abs(maxPad - (*mp2)->GetPad()) < fMaxPadRange && //4\r
 792           \r
 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
 796           \r
 797           clusterValue+=(mp2)->GetSum();\r
 798 \r
 799           nUsedPads++;\r
 800           nUsedRows.push_back((mp2)->GetRow());\r
 801 \r
 802           AliSimDigits *digrowTmp1;\r
 803           if(fRawData){\r
 804             digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());\r
 805           }else{\r
 806             digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());\r
 807           }\r
 808           digrowTmp1->ExpandBuffer(); //decrunch\r
 809           \r
 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
 815           }\r
 816           \r
 817           tmp->SetCenter();//set center of the cluster\r
 818 \r
 819           //which one is bigger\r
 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
 837             }\r
 838           }\r
 839           maximaInSector->RemoveAt(it2);\r
 840           maximaInSector->Compress();\r
 841           it2--;\r
 842 \r
 843           //maximaInSector.erase(mp2);\r
 844           //mp2--;\r
 845         }\r
 846       }//inside loop\r
 847 \r
 848       tmp->SetSize();\r
 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
 854       //if(clusterValue/tmp->fCluster.size()<fValueToSize)continue;\r
 855       if(clusterValue/(tmp->GetSize())<fValueToSize)continue;\r
 856 \r
 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
 870 \r
 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
 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
 897 void AliTPCclustererKr::GetXY(Short_t sec,Short_t row,Short_t pad,Double_t& xGlob,Double_t& yGlob){\r
 898   //\r
 899   //gives global XY coordinate of the pad\r
 900   //\r
 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
 908   Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad\r
 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
 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
 926 \r
 927   return;\r
 928 }\r