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