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 cluster class
  20 //
  21 // Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl
  22 //-----------------------------------------------------------------
  23
  24 #include "AliTPCclustererKr.h"
  25 #include "AliTPCclusterKr.h"
  26 #include <vector>
  27 #include "TObject.h"
  28 #include "AliPadMax.h"
  29 #include "AliSimDigits.h"
  30 #include "AliTPCv4.h"
  31 #include "AliTPCParam.h"
  32 #include "AliTPCDigitsArray.h"
  33 #include "AliTPCvtpr.h"
  34 #include "AliTPCClustersRow.h"
  35 #include "TTree.h"
  36
  37 //used in raw data finder
  38 #include "AliTPCROC.h"
  39 #include "AliTPCCalPad.h"
  40 #include "AliTPCAltroMapping.h"
  41 #include "AliTPCcalibDB.h"
  42 #include "AliTPCRawStream.h"
  43 #include "AliTPCRecoParam.h"
  44 #include "AliTPCReconstructor.h"
  45 #include "AliRawReader.h"
  46 #include "AliTPCCalROC.h"
  47
  48 ClassImp(AliTPCclustererKr)
  49
  50
  51 AliTPCclustererKr::AliTPCclustererKr()
  52   :TObject(),
  53   fRawData(kFALSE),
  54   fRowCl(0),
  55   fInput(0),
  56   fOutput(0),
  57   fParam(0),
  58   fDigarr(0),
  59   fRecoParam(0),
  60   fIsOldRCUFormat(kFALSE)
  61 {
  62 //
  63 // default constructor
  64 //
  65 }
  66
  67 AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr &param)
  68   :TObject(),
  69   fRawData(kFALSE),
  70   fRowCl(0),
  71   fInput(0),
  72   fOutput(0),
  73   fParam(0),
  74   fDigarr(0),
  75   fRecoParam(0),
  76   fIsOldRCUFormat(kFALSE)
  77 {
  78 //
  79 // copy constructor
  80 //
  81   fParam = param.fParam;
  82   fRecoParam=param.fRecoParam;
  83   fIsOldRCUFormat=param.fIsOldRCUFormat;
  84   fRawData=param.fRawData;
  85   fRowCl =param.fRowCl ;
  86   fInput =param.fInput ;
  87   fOutput=param.fOutput;
  88   fDigarr=param.fDigarr;
  89 }
  90
  91 AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
  92 {
  93   fParam = param.fParam;
  94   fRecoParam=param.fRecoParam;
  95   fIsOldRCUFormat=param.fIsOldRCUFormat;
  96   fRawData=param.fRawData;
  97   fRowCl =param.fRowCl ;
  98   fInput =param.fInput ;
  99   fOutput=param.fOutput;
 100   fDigarr=param.fDigarr;
 101   return (*this);
 102 }
 103
 104 AliTPCclustererKr::~AliTPCclustererKr()
 105 {
 106   //
 107   // destructor
 108   //
 109 }
 110
 111 void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
 112 {
 113   if (recoParam) {
 114     fRecoParam = recoParam;
 115   }else{
 116     //set default parameters if not specified
 117     fRecoParam = AliTPCReconstructor::GetRecoParam();
 118     if (!fRecoParam)  fRecoParam = AliTPCRecoParam::GetLowFluxParam();
 119   }
 120   return;
 121 }
 122
 123
 124 ////____________________________________________________________________________
 125 ////       I/O
 126 void AliTPCclustererKr::SetInput(TTree * tree)
 127 {
 128   //
 129   // set input tree with digits
 130   //
 131   fInput = tree;
 132   if  (!fInput->GetBranch("Segment")){
 133     cerr<<"AliTPCclusterKr::findClusterKr(): no proper input tree !\n";
 134     fInput=0;
 135     return;
 136   }
 137 }
 138
 139 void AliTPCclustererKr::SetOutput(TTree * tree)
 140 {
 141   //
 142   // set output
 143   //
 144   fOutput= tree;
 145   AliTPCClustersRow clrow;
 146   AliTPCClustersRow *pclrow=&clrow;
 147   clrow.SetClass("AliTPCclusterKr");
 148   clrow.SetArray(1); // to make Clones array
 149   fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
 150 }
 151
 152 ////____________________________________________________________________________
 153 //// with new I/O
 154 Int_t AliTPCclustererKr::finderIO()
 155 {
 156   // Krypton cluster finder for simulated events from MC
 157
 158   if (!fInput) {
 159     Error("Digits2Clusters", "input tree not initialised");
 160     return 10;
 161   }
 162
 163   if (!fOutput) {
 164     Error("Digits2Clusters", "output tree not initialised");
 165     return 11;
 166   }
 167
 168   findClusterKrIO();
 169   return 0;
 170 }
 171
 172 Int_t AliTPCclustererKr::finderIO(AliRawReader* rawReader)
 173 {
 174   // Krypton cluster finder for the TPC raw data
 175   //
 176   // fParam must be defined before
 177
 178   // consider noiceROC or not
 179
 180   if(rawReader)fRawData=kTRUE; //set flag to data
 181
 182   if (!fOutput) {
 183     Error("Digits2Clusters", "output tree not initialised");
 184     return 11;
 185   }
 186
 187   fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
 188   //   used later for memory allocation
 189
 190   Bool_t isAltro=kFALSE;
 191
 192   AliTPCROC * roc = AliTPCROC::Instance();
 193   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 194   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
 195   //
 196   AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
 197
 198   const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
 199   const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
 200   const Int_t kNS = kNIS + kNOS;//all sectors
 201   Int_t zeroSup = fParam->GetZeroSup();//zero suppression parameter
 202
 203   //crate TPC view
 204   AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
 205   digarr->Setup(fParam);//as usually parameters
 206
 207   //
 208   // Loop over sectors
 209   //
 210   for(Int_t sec = 0; sec < kNS; sec++) {
 211     AliTPCCalROC * noiseROC;
 212     if(noiseTPC==0x0){
 213       noiseROC =new AliTPCCalROC(sec);//noise=0
 214     }
 215     else{
 216       noiseROC = noiseTPC->GetCalROC(sec);  // noise per given sector
 217     }
 218     Int_t nRows = 0; //number of rows in sector
 219     Int_t nDDLs = 0; //number of DDLs
 220     Int_t indexDDL = 0; //DDL index
 221     if (sec < kNIS) {
 222       nRows = fParam->GetNRowLow();
 223       nDDLs = 2;
 224       indexDDL = sec * 2;
 225     }else {
 226       nRows = fParam->GetNRowUp();
 227       nDDLs = 4;
 228       indexDDL = (sec-kNIS) * 4 + kNIS * 2;
 229     }
 230
 231     //
 232     // Load the raw data for corresponding DDLs
 233     //
 234     rawReader->Reset();
 235     input.SetOldRCUFormat(fIsOldRCUFormat);
 236     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
 237
 238     if(input.Next()) {
 239       isAltro=kTRUE;
 240       // Allocate memory for rows in sector (pads(depends on row) x timebins)
 241       for(Int_t row = 0; row < nRows; row++) {
 242         digarr->CreateRow(sec,row);
 243       }//end loop over rows
 244     }
 245     rawReader->Reset();
 246     input.SetOldRCUFormat(fIsOldRCUFormat);
 247     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
 248
 249     //
 250     // Begin loop over altro data
 251     //
 252     while (input.Next()) {
 253
 254       //check sector consistency
 255       if (input.GetSector() != sec)
 256         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",sec,input.GetSector()));
 257
 258       //check row consistency
 259       Short_t iRow = input.GetRow();
 260       if (iRow < 0 || iRow >= nRows){
 261         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
 262                       iRow, 0, nRows -1));
 263         continue;
 264       }
 265
 266       //check pad consistency
 267       Short_t iPad = input.GetPad();
 268       if (iPad < 0 || iPad >= (Short_t)(roc->GetNPads(sec,iRow))) {
 269         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
 270                       iPad, 0, roc->GetNPads(sec,iRow) ));
 271         continue;
 272       }
 273
 274       //check time consistency
 275       Short_t iTimeBin = input.GetTime();
 276       if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
 277         //cout<<iTimeBin<<endl;
 278         continue;
 279         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
 280                       iTimeBin, 0, fRecoParam->GetLastBin() -1));
 281       }
 282
 283       //signal
 284       Int_t signal = input.GetSignal();
 285       Short_t timeWindowLow=60, timeWindowHigh=950;
 286       if (signal <= zeroSup ||
 287           iTimeBin < timeWindowLow ||
 288           iTimeBin > timeWindowHigh
 289           ) {
 290         digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
 291         continue;
 292       }
 293       if (!noiseROC) continue;
 294       Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
 295       if (noiseOnPad>2) continue; // consider noisy pad as dead
 296
 297       if(signal<4*noiseOnPad){
 298         digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
 299         continue;
 300       }
 301
 302 //      if(signal>12){
 303 //      cout<<" s "<<sec<<" r "<<iRow<<" p "<<iPad<<" t "<<iTimeBin<<" dig "<<signal<<endl;
 304       digarr->GetRow(sec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
 305 //      }else digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
 306
 307
 308     }//end of loop over altro data
 309   }//end of loop over sectors
 310
 311   SetDigArr(digarr);
 312   if(isAltro) findClusterKrIO();
 313   delete digarr;
 314
 315   return 0;
 316 }
 317
 318 ////____________________________________________________________________________
 319 Int_t AliTPCclustererKr::findClusterKrIO()
 320 {
 321   //fParam and  fDigarr must be set to run this method
 322
 323   Int_t cluster_counter=0;
 324   const Short_t NTotalSector=fParam->GetNSector();//number of sectors
 325   for(Short_t sec=0; sec<NTotalSector; sec++){
 326
 327     //vector of maxima for each sector
 328     std::vector<AliPadMax*> maxima_in_sector;
 329
 330     const Short_t Nrows=fParam->GetNRow(sec);//number of rows in sector
 331     for(Short_t row=0; row<Nrows; row++){
 332       AliSimDigits *digrow;
 333       if(fRawData){
 334         digrow = (AliSimDigits*)fDigarr->GetRow(sec,row);//real data
 335       }else{
 336         digrow = (AliSimDigits*)fDigarr->LoadRow(sec,row);//MC
 337       }
 338       if(digrow){//if pointer exist
 339         digrow->ExpandBuffer(); //decrunch
 340         const Short_t npads = digrow->GetNCols();  // number of pads
 341         const Short_t ntime = digrow->GetNRows(); // number of timebins
 342         for(Short_t np=0;np<npads;np++){
 343
 344           Short_t tb_max=-1;//timebin of maximum
 345           Short_t value_maximum=-1;//value of maximum in adc
 346           Short_t increase_begin=-1;//timebin when increase starts
 347           Short_t sum_adc=0;//sum of adc on the pad in maximum surrounding
 348           bool if_increase_begin=true;//flag - check if increasing started
 349           bool if_maximum=false;//flag - check if it could be maximum
 350
 351           for(Short_t nt=0;nt<ntime;nt++){
 352             Short_t adc = digrow->GetDigitFast(nt,np);
 353             if(adc<6){//standard was 3
 354               if(if_maximum){
 355                 if(nt-1-increase_begin<1){//at least 2 time bins, was = 1
 356                   tb_max=-1;
 357                   value_maximum=-1;
 358                   increase_begin=-1;
 359                   sum_adc=0;
 360                   if_increase_begin=true;
 361                   if_maximum=false;
 362                   continue;
 363                 }
 364                 //insert maximum, default values and set flags
 365                 AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
 366                                                                   tb_max,
 367                                                                   np,
 368                                                                   row),
 369                                                        increase_begin,
 370                                                        nt-1,
 371                                                        sum_adc);
 372                 maxima_in_sector.push_back(one_maximum);
 373
 374                 tb_max=-1;
 375                 value_maximum=-1;
 376                 increase_begin=-1;
 377                 sum_adc=0;
 378                 if_increase_begin=true;
 379                 if_maximum=false;
 380               }
 381               continue;
 382             }
 383
 384             if(if_increase_begin){
 385               if_increase_begin=false;
 386               increase_begin=nt;
 387             }
 388
 389             if(adc>value_maximum){
 390               tb_max=nt;
 391               value_maximum=adc;
 392               if_maximum=true;
 393             }
 394             sum_adc+=adc;
 395             if(nt==ntime-1 && if_maximum && ntime-1-increase_begin>2){//on the edge
 396               //insert maximum, default values and set flags
 397               AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
 398                                                                 tb_max,
 399                                                                 np,
 400                                                                 row),
 401                                                      increase_begin,
 402                                                      nt-1,
 403                                                      sum_adc);
 404               maxima_in_sector.push_back(one_maximum);
 405
 406               tb_max=-1;
 407               value_maximum=-1;
 408               increase_begin=-1;
 409               sum_adc=0;
 410               if_increase_begin=true;
 411               if_maximum=false;
 412               continue;
 413             }
 414
 415           }//end loop over timebins
 416         }//end loop over pads
 417 //      }else{
 418 //      cout<<"Pointer does not exist!!"<<endl;
 419       }//end if poiner exists
 420     }//end loop over rows
 421
 422
 423     Short_t max_dig=0;
 424     Short_t max_sum_adc=0;
 425     Short_t max_nt=0;
 426     Short_t max_np=0;
 427     Short_t max_row=0;
 428
 429     for( std::vector<AliPadMax*>::iterator mp1  = maxima_in_sector.begin();
 430          mp1 != maxima_in_sector.end(); ++mp1 ) {
 431
 432       AliTPCclusterKr *tmp=new AliTPCclusterKr();
 433
 434       Short_t n_used_pads=1;
 435       Int_t cluster_value=0;
 436       cluster_value+=(*mp1)->GetSum();
 437
 438       max_dig      =(*mp1)->GetAdc() ;
 439       max_sum_adc  =(*mp1)->GetSum() ;
 440       max_nt       =(*mp1)->GetTime();
 441       max_np       =(*mp1)->GetPad() ;
 442       max_row      =(*mp1)->GetRow() ;
 443
 444       AliSimDigits *digrow_tmp;
 445       if(fRawData){
 446         digrow_tmp = (AliSimDigits*)fDigarr->GetRow(sec,(*mp1)->GetRow());
 447       }else{
 448         digrow_tmp = (AliSimDigits*)fDigarr->LoadRow(sec,(*mp1)->GetRow());
 449       }
 450
 451       digrow_tmp->ExpandBuffer(); //decrunch
 452       for(Short_t tb=(*mp1)->GetBegin(); tb<((*mp1)->GetEnd())+1; tb++){
 453         Short_t adc_tmp = digrow_tmp->GetDigitFast(tb,(*mp1)->GetPad());
 454         AliTPCvtpr *vtpr=new AliTPCvtpr(adc_tmp,tb,(*mp1)->GetPad(),(*mp1)->GetRow());
 455         tmp->fCluster.push_back(vtpr);
 456       }
 457
 458       maxima_in_sector.erase(mp1);
 459       mp1--;
 460
 461       for( std::vector<AliPadMax*>::iterator mp2  = maxima_in_sector.begin();
 462            mp2 != maxima_in_sector.end(); ++mp2 ) {
 463
 464         Short_t two_pad_max=4;
 465         Short_t two_row_max=3;
 466
 467         if(abs(max_row - (*mp2)->GetRow())<two_row_max &&
 468            abs(max_np - (*mp2)->GetPad())<two_pad_max  &&
 469            abs(max_nt - (*mp2)->GetTime())<7){//was 7
 470
 471           cluster_value+=(*mp2)->GetSum();
 472
 473           n_used_pads++;
 474
 475           AliSimDigits *digrow_tmp1;
 476           if(fRawData){
 477             digrow_tmp1 = (AliSimDigits*)fDigarr->GetRow(sec,(*mp2)->GetRow());
 478           }else{
 479             digrow_tmp1 = (AliSimDigits*)fDigarr->LoadRow(sec,(*mp2)->GetRow());
 480           }
 481           digrow_tmp1->ExpandBuffer(); //decrunch
 482
 483           for(Short_t tb=(*mp2)->GetBegin(); tb<(*mp2)->GetEnd()+1; tb++){
 484             Short_t adc_tmp = digrow_tmp1->GetDigitFast(tb,(*mp2)->GetPad());
 485             AliTPCvtpr *vtpr=new AliTPCvtpr(adc_tmp,tb,(*mp2)->GetPad(),(*mp2)->GetRow());
 486               tmp->fCluster.push_back(vtpr);
 487           }
 488
 489           //which one is bigger
 490           if( (*mp2)->GetAdc() > max_dig ){
 491             max_dig      =(*mp2)->GetAdc() ;
 492             max_sum_adc  =(*mp2)->GetSum() ;
 493             max_nt       =(*mp2)->GetTime();
 494             max_np       =(*mp2)->GetPad() ;
 495             max_row      =(*mp2)->GetRow() ;
 496           } else if ( (*mp2)->GetAdc() == max_dig ){
 497             if( (*mp2)->GetSum() > max_sum_adc){
 498               max_dig      =(*mp2)->GetAdc() ;
 499               max_sum_adc  =(*mp2)->GetSum() ;
 500               max_nt       =(*mp2)->GetTime();
 501               max_np       =(*mp2)->GetPad() ;
 502               max_row      =(*mp2)->GetRow() ;
 503             }
 504           }
 505           maxima_in_sector.erase(mp2);
 506           mp2--;
 507         }
 508       }//inside loop
 509
 510       //through out ADC=6,7 on 1 pad, 2 tb and ADC=12 on 2 pads,2 tb
 511       if(n_used_pads==1 && cluster_value/tmp->fCluster.size()<3.6)continue;
 512       if(n_used_pads==2 && cluster_value/tmp->fCluster.size()<3.1)continue;
 513
 514       tmp->SetADCcluster(cluster_value);
 515       tmp->SetNpads(n_used_pads);
 516       tmp->SetMax(AliTPCvtpr(max_dig,max_nt,max_np,max_row));
 517       tmp->SetSec(sec);
 518       tmp->SetSize();
 519
 520       cluster_counter++;
 521
 522       //save each cluster into file
 523
 524       AliTPCClustersRow *clrow= new AliTPCClustersRow();
 525       fRowCl=clrow;
 526       clrow->SetClass("AliTPCclusterKr");
 527       clrow->SetArray(1);
 528       fOutput->GetBranch("Segment")->SetAddress(&clrow);
 529
 530       Int_t tmp_cluster=0;
 531       TClonesArray * arr = fRowCl->GetArray();
 532       AliTPCclusterKr * cl = new ((*arr)[tmp_cluster]) AliTPCclusterKr(*tmp);
 533
 534       fOutput->Fill();
 535       delete clrow;
 536
 537       //end of save each cluster into file adc.root
 538     }//outer loop
 539
 540   }//end sector for
 541   cout<<"Number of clusters in event: "<<cluster_counter<<endl;
 542
 543   return 0;
 544 }
 545
 546 ////____________________________________________________________________________