TPC/AliTPCclustererKr.cxx

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