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