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   return 0;
 624 }
 625 ////____________________________________________________________________________
 626 Int_t AliTPCclustererKr::FindClusterKrIO()
 627 {
 628
 629   //
 630   //fParam and  fDigarr must be set to run this method
 631   //
 632
 633   Int_t clusterCounter=0;
 634   const Int_t nTotalSector=fParam->GetNSector();//number of sectors
 635   for(Int_t iSec=0; iSec<nTotalSector; ++iSec){
 636     CleanSector(iSec);
 637     //vector of maxima for each sector
 638     //std::vector<AliPadMax*> maximaInSector;
 639     TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*
 640
 641     //
 642     //  looking for the maxima on the pad
 643     //
 644
 645
 646
 647     const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector
 648     for(Int_t iRow=0; iRow<kNRows; ++iRow){
 649       AliSimDigits *digrow;
 650       if(fRawData){
 651         digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
 652       }else{
 653         digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
 654       }
 655       if(digrow){//if pointer exist
 656         digrow->ExpandBuffer(); //decrunch
 657         const Int_t kNPads = digrow->GetNCols();  // number of pads
 658         const Int_t kNTime = digrow->GetNRows(); // number of timebins
 659         for(Int_t iPad=0;iPad<kNPads;iPad++){
 660
 661           Int_t timeBinMax=-1;//timebin of maximum
 662           Int_t valueMaximum=-1;//value of maximum in adc
 663           Int_t increaseBegin=-1;//timebin when increase starts
 664           Int_t sumAdc=0;//sum of adc on the pad in maximum surrounding
 665           bool ifIncreaseBegin=true;//flag - check if increasing started
 666           bool ifMaximum=false;//flag - check if it could be maximum
 667           Short_t* val = digrow->GetDigitsColumn(iPad);
 668           for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
 669             if (!ifMaximum)  {
 670               if (val[iTimeBin]==-1) break;   // 0 until the end
 671               for( ; iTimeBin<kNTime-2&&val[iTimeBin]<fMinAdc ;iTimeBin++){}
 672             }
 673             //
 674             Short_t adc = val[iTimeBin];
 675             if(adc<fMinAdc){//standard was 3
 676               if(ifMaximum){
 677                 if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
 678                   timeBinMax=-1;
 679                   valueMaximum=-1;
 680                   increaseBegin=-1;
 681                   sumAdc=0;
 682                   ifIncreaseBegin=true;
 683                   ifMaximum=false;
 684                   continue;
 685                 }
 686                 //insert maximum, default values and set flags
 687                 //Double_t xCord,yCord;
 688                 //GetXY(iSec,iRow,iPad,xCord,yCord);
 689                 Double_t x[]={iRow,iPad,iTimeBin};
 690                 Int_t i[]={iSec};
 691                 AliTPCTransform trafo;
 692                 trafo.Transform(x,i,0,1);
 693
 694                 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
 695                                                                  timeBinMax,
 696                                                                  iPad,
 697                                                                  iRow,
 698                                                                  x[0],//xCord,
 699                                                                  x[1],//yCord,
 700                                                                  x[2]/*timeBinMax*/),
 701                                                       increaseBegin,
 702                                                       iTimeBin-1,
 703                                                       sumAdc);
 704                 maximaInSector->AddLast(oneMaximum);
 705
 706                 timeBinMax=-1;
 707                 valueMaximum=-1;
 708                 increaseBegin=-1;
 709                 sumAdc=0;
 710                 ifIncreaseBegin=true;
 711                 ifMaximum=false;
 712               }
 713               continue;
 714             }
 715
 716             if(ifIncreaseBegin){
 717               ifIncreaseBegin=false;
 718               increaseBegin=iTimeBin;
 719             }
 720
 721             if(adc>valueMaximum){
 722               timeBinMax=iTimeBin;
 723               valueMaximum=adc;
 724               ifMaximum=true;
 725             }
 726             sumAdc+=adc;
 727             if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
 728               //at least 3 timebins
 729               //insert maximum, default values and set flags
 730               //Double_t xCord,yCord;
 731               //GetXY(iSec,iRow,iPad,xCord,yCord);
 732               Double_t x[]={iRow,iPad,iTimeBin};
 733               Int_t i[]={iSec};
 734               AliTPCTransform trafo;
 735               trafo.Transform(x,i,0,1);
 736               AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
 737                                                                timeBinMax,
 738                                                                iPad,
 739                                                                iRow,
 740                                                                x[0],//xCord,
 741                                                                x[1],//yCord,
 742                                                                x[2]/*timeBinMax*/),
 743                                                     increaseBegin,
 744                                                     iTimeBin-1,
 745                                                     sumAdc);
 746               maximaInSector->AddLast(oneMaximum);
 747
 748               timeBinMax=-1;
 749               valueMaximum=-1;
 750               increaseBegin=-1;
 751               sumAdc=0;
 752               ifIncreaseBegin=true;
 753               ifMaximum=false;
 754               continue;
 755             }
 756
 757           }//end loop over timebins
 758         }//end loop over pads
 759 //      }else{
 760 //      cout<<"Pointer does not exist!!"<<endl;
 761       }//end if poiner exists
 762     }//end loop over rows
 763
 764     MakeClusters(maximaInSector,iSec,clusterCounter);
 765     //
 766     maximaInSector->SetOwner(kTRUE);
 767     maximaInSector->Delete();
 768     delete maximaInSector;
 769   }//end sector for
 770   cout<<"Number of clusters in event: "<<clusterCounter<<endl;
 771   return 0;
 772 }
 773
 774 void AliTPCclustererKr::MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter){
 775   //
 776   // Make clusters
 777   //
 778
 779   Int_t maxDig=0;
 780   Int_t maxSumAdc=0;
 781   Int_t maxTimeBin=0;
 782   Int_t maxPad=0;
 783   Int_t maxRow=0;
 784   Double_t maxX=0;
 785   Double_t maxY=0;
 786   Double_t maxT=0;
 787   Int_t entriesArr = maximaInSector->GetEntriesFast();
 788   for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) {
 789
 790     AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1);
 791     if (!mp1) continue;
 792     AliTPCclusterKr clusterKr;
 793
 794     Int_t nUsedPads=1;
 795     Int_t clusterValue=0;
 796     clusterValue+=(mp1)->GetSum();
 797     list<Int_t> nUsedRows;
 798     nUsedRows.push_back((mp1)->GetRow());
 799
 800     maxDig      =(mp1)->GetAdc() ;
 801     maxSumAdc   =(mp1)->GetSum() ;
 802     maxTimeBin  =(mp1)->GetTime();
 803     maxPad      =(mp1)->GetPad() ;
 804     maxRow      =(mp1)->GetRow() ;
 805     maxX        =(mp1)->GetX();
 806     maxY        =(mp1)->GetY();
 807     maxT        =(mp1)->GetT();
 808
 809     AliSimDigits *digrowTmp;
 810     if(fRawData){
 811       digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());
 812     }else{
 813       digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());
 814     }
 815
 816     digrowTmp->ExpandBuffer(); //decrunch
 817
 818     for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){
 819       Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad());
 820       AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT());
 821       clusterKr.AddDigitToCluster(vtpr);
 822     }
 823     clusterKr.SetCenter();//set centr of the cluster
 824
 825     for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) {
 826       AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2);
 827       if (!mp2) continue;
 828       if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX())>fMaxPadRangeCm) continue;
 829       if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue;
 830       if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) >fMaxTimeRange) continue;
 831
 832       {
 833         clusterValue+=(mp2)->GetSum();
 834
 835         nUsedPads++;
 836         nUsedRows.push_back((mp2)->GetRow());
 837
 838         AliSimDigits *digrowTmp1;
 839         if(fRawData){
 840           digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());
 841         }else{
 842           digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());
 843         }
 844         digrowTmp1->ExpandBuffer(); //decrunch
 845
 846         for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){
 847           Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad());
 848           AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT());
 849           clusterKr.AddDigitToCluster(vtpr);
 850         }
 851
 852         clusterKr.SetCenter();//set center of the cluster
 853
 854         //which one is bigger
 855         if( (mp2)->GetAdc() > maxDig ){
 856           maxDig      =(mp2)->GetAdc() ;
 857           maxSumAdc   =(mp2)->GetSum() ;
 858           maxTimeBin  =(mp2)->GetTime();
 859           maxPad      =(mp2)->GetPad() ;
 860           maxRow      =(mp2)->GetRow() ;
 861           maxX        =(mp2)->GetX() ;
 862           maxY        =(mp2)->GetY() ;
 863           maxT        =(mp2)->GetT() ;
 864         } else if ( (mp2)->GetAdc() == maxDig ){
 865           if( (mp2)->GetSum() > maxSumAdc){
 866             maxDig      =(mp2)->GetAdc() ;
 867             maxSumAdc   =(mp2)->GetSum() ;
 868             maxTimeBin  =(mp2)->GetTime();
 869             maxPad      =(mp2)->GetPad() ;
 870             maxRow      =(mp2)->GetRow() ;
 871             maxX        =(mp2)->GetX() ;
 872             maxY        =(mp2)->GetY() ;
 873             maxT        =(mp2)->GetT() ;
 874           }
 875         }
 876         delete maximaInSector->RemoveAt(it2);
 877       }
 878     }//inside loop
 879     delete maximaInSector->RemoveAt(it1);
 880     clusterKr.SetSize();
 881     //through out clusters on the edge and noise
 882     //if(clusterValue/clusterKr.fCluster.size()<fValueToSize)continue;
 883     if(clusterValue/(clusterKr.GetSize())<fValueToSize)continue;
 884
 885     clusterKr.SetADCcluster(clusterValue);
 886     clusterKr.SetNPads(nUsedPads);
 887     clusterKr.SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxT));
 888     clusterKr.SetSec(iSec);
 889     clusterKr.SetSize();
 890
 891     nUsedRows.sort();
 892     nUsedRows.unique();
 893     clusterKr.SetNRows(nUsedRows.size());
 894     clusterKr.SetCenter();
 895
 896     clusterCounter++;
 897
 898
 899     //save each cluster into file
 900     if (fOutput){
 901       (*fOutput)<<"Kr"<<
 902         "Cl.="<<&clusterKr<<
 903         "\n";
 904     }
 905     //end of save each cluster into file adc.root
 906   }//outer loop
 907 }
 908
 909
 910
 911 ////____________________________________________________________________________
 912
 913
 914 void AliTPCclustererKr::GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob){
 915   //
 916   //gives global XY coordinate of the pad
 917   //
 918
 919   Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
 920
 921   Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
 922   Float_t padXSize;
 923   if(sec<fParam->GetNInnerSector())padXSize=0.4;
 924   else padXSize=0.6;
 925   Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad
 926
 927   Float_t sin,cos;
 928   fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
 929
 930   Double_t xGlob1 =  xLocal * cos + yLocal * sin;
 931   Double_t yGlob1 = -xLocal * sin + yLocal * cos;
 932
 933
 934   Double_t rot=0;
 935   rot=TMath::Pi()/2.;
 936
 937   xGlob =  xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);
 938   yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);
 939
 940    yGlob=-1*yGlob;
 941    if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;
 942
 943
 944   return;
 945 }