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