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