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   = noiseTPC->GetCalROC(sec);  // noise per given sector
 212
 213     Int_t nRows = 0; //number of rows in sector
 214     Int_t nDDLs = 0; //number of DDLs
 215     Int_t indexDDL = 0; //DDL index
 216     if (sec < kNIS) {
 217       nRows = fParam->GetNRowLow();
 218       nDDLs = 2;
 219       indexDDL = sec * 2;
 220     }else {
 221       nRows = fParam->GetNRowUp();
 222       nDDLs = 4;
 223       indexDDL = (sec-kNIS) * 4 + kNIS * 2;
 224     }
 225
 226     //
 227     // Load the raw data for corresponding DDLs
 228     //
 229     rawReader->Reset();
 230     input.SetOldRCUFormat(fIsOldRCUFormat);
 231     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
 232
 233     if(input.Next()) {
 234       isAltro=kTRUE;
 235       // Allocate memory for rows in sector (pads(depends on row) x timebins)
 236       for(Int_t row = 0; row < nRows; row++) {
 237         digarr->CreateRow(sec,row);
 238       }//end loop over rows
 239     }
 240     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
 241
 242     //
 243     // Begin loop over altro data
 244     //
 245     while (input.Next()) {
 246
 247       //check sector consistency
 248       if (input.GetSector() != sec)
 249         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",sec,input.GetSector()));
 250
 251       //check row consistency
 252       Short_t iRow = input.GetRow();
 253       if (iRow < 0 || iRow >= nRows){
 254         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
 255                       iRow, 0, nRows -1));
 256         continue;
 257       }
 258
 259       //check pad consistency
 260       Short_t iPad = input.GetPad();
 261       if (iPad < 0 || iPad >= (Short_t)(roc->GetNPads(sec,iRow))) {
 262         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
 263                       iPad, 0, roc->GetNPads(sec,iRow) ));
 264         continue;
 265       }
 266
 267       //check time consistency
 268       Short_t iTimeBin = input.GetTime();
 269       if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
 270         //cout<<iTimeBin<<endl;
 271         continue;
 272         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
 273                       iTimeBin, 0, fRecoParam->GetLastBin() -1));
 274       }
 275
 276       //signal
 277       Int_t signal = input.GetSignal();
 278       if (signal <= zeroSup) {
 279         //continue;
 280         digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
 281       }
 282       (//(AliSimDigits*)
 283        digarr->GetRow(sec,iRow))->SetDigitFast(signal,iTimeBin,iPad);
 284     }//end of loop over altro data
 285   }//end of loop over sectors
 286
 287   SetDigArr(digarr);
 288   if(isAltro) findClusterKrIO();
 289   delete digarr;
 290
 291   return 0;
 292 }
 293
 294 ////____________________________________________________________________________
 295 Int_t AliTPCclustererKr::findClusterKrIO()
 296 {
 297   //fParam and  fDigarr must be set to run this method
 298
 299   Int_t cluster_counter=0;
 300   const Short_t NTotalSector=fParam->GetNSector();//number of sectors
 301   for(Short_t sec=0; sec<NTotalSector; sec++){
 302
 303     //vector of maxima for each sector
 304     std::vector<AliPadMax*> maxima_in_sector;
 305
 306     const Short_t Nrows=fParam->GetNRow(sec);//number of rows in sector
 307     for(Short_t row=0; row<Nrows; row++){
 308       AliSimDigits *digrow;
 309       if(fRawData){
 310         digrow = (AliSimDigits*)fDigarr->GetRow(sec,row);//real data
 311       }else{
 312         digrow = (AliSimDigits*)fDigarr->LoadRow(sec,row);//MC
 313       }
 314       if(digrow){//if pointer exist
 315         digrow->ExpandBuffer(); //decrunch
 316         const Short_t npads = digrow->GetNCols();  // number of pads
 317         const Short_t ntime = digrow->GetNRows(); // number of timebins
 318         for(Short_t np=0;np<npads;np++){
 319
 320           Short_t tb_max=-1;//timebin of maximum
 321           Short_t value_maximum=-1;//value of maximum in adc
 322           Short_t increase_begin=-1;//timebin when increase starts
 323           Short_t sum_adc=0;//sum of adc on the pad in maximum surrounding
 324           bool if_increase_begin=true;//flag - check if increasing start
 325           bool if_maximum=false;//flag - check if it could be maximum
 326
 327           for(Short_t nt=0;nt<ntime;nt++){
 328             Short_t adc = digrow->GetDigitFast(nt,np);
 329             if(adc<3){
 330               if(if_maximum){
 331                 if(nt-1-increase_begin<1){//at least 2 time bins
 332                   tb_max=-1;
 333                   value_maximum=-1;
 334                   increase_begin=-1;
 335                   sum_adc=0;
 336                   if_increase_begin=true;
 337                   if_maximum=false;
 338                   continue;
 339                 }
 340                 //insert maximum, default values and set flags
 341                 AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
 342                                                                   tb_max,
 343                                                                   np,
 344                                                                   row),
 345                                                        increase_begin,
 346                                                        nt-1,
 347                                                        sum_adc);
 348                 maxima_in_sector.push_back(one_maximum);
 349
 350                 tb_max=-1;
 351                 value_maximum=-1;
 352                 increase_begin=-1;
 353                 sum_adc=0;
 354                 if_increase_begin=true;
 355                 if_maximum=false;
 356               }
 357               continue;
 358             }
 359
 360             if(if_increase_begin){
 361               if_increase_begin=false;
 362               increase_begin=nt;
 363             }
 364
 365             if(adc>value_maximum){
 366               tb_max=nt;
 367               value_maximum=adc;
 368               if_maximum=true;
 369             }
 370             sum_adc+=adc;
 371             if(nt==ntime-1 && if_maximum && ntime-1-increase_begin>2){//on the edge
 372               //insert maximum, default values and set flags
 373               AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
 374                                                                 tb_max,
 375                                                                 np,
 376                                                                 row),
 377                                                      increase_begin,
 378                                                      nt-1,
 379                                                      sum_adc);
 380               maxima_in_sector.push_back(one_maximum);
 381
 382               tb_max=-1;
 383               value_maximum=-1;
 384               increase_begin=-1;
 385               sum_adc=0;
 386               if_increase_begin=true;
 387               if_maximum=false;
 388               continue;
 389             }
 390
 391           }//end loop over timebins
 392         }//end loop over pads
 393 //      }else{
 394 //      cout<<"Pointer does not exist!!"<<endl;
 395       }//end if poiner exists
 396     }//end loop over rows
 397
 398
 399     Short_t max_dig=0;
 400     Short_t max_sum_adc=0;
 401     Short_t max_nt=0;
 402     Short_t max_np=0;
 403     Short_t max_row=0;
 404
 405     for( std::vector<AliPadMax*>::iterator mp1  = maxima_in_sector.begin();
 406          mp1 != maxima_in_sector.end(); ++mp1 ) {
 407
 408       AliTPCclusterKr *tmp=new AliTPCclusterKr();
 409
 410       Short_t n_used_pads=1;
 411       Short_t cluster_value=0;
 412       cluster_value+=(*mp1)->GetSum();
 413
 414       max_dig      =(*mp1)->GetAdc() ;
 415       max_sum_adc  =(*mp1)->GetSum() ;
 416       max_nt       =(*mp1)->GetTime();
 417       max_np       =(*mp1)->GetPad() ;
 418       max_row      =(*mp1)->GetRow() ;
 419
 420       AliSimDigits *digrow_tmp;
 421       if(fRawData){
 422         digrow_tmp = (AliSimDigits*)fDigarr->GetRow(sec,(*mp1)->GetRow());
 423       }else{
 424         digrow_tmp = (AliSimDigits*)fDigarr->LoadRow(sec,(*mp1)->GetRow());
 425       }
 426
 427       digrow_tmp->ExpandBuffer(); //decrunch
 428       for(Short_t tb=(*mp1)->GetBegin(); tb<((*mp1)->GetEnd())+1; tb++){
 429         Short_t adc_tmp = digrow_tmp->GetDigitFast(tb,(*mp1)->GetPad());
 430         AliTPCvtpr *vtpr=new AliTPCvtpr(adc_tmp,tb,(*mp1)->GetPad(),(*mp1)->GetRow());
 431         tmp->fCluster.push_back(vtpr);
 432       }
 433
 434       maxima_in_sector.erase(mp1);
 435       mp1--;
 436
 437       for( std::vector<AliPadMax*>::iterator mp2  = maxima_in_sector.begin();
 438            mp2 != maxima_in_sector.end(); ++mp2 ) {
 439
 440         Short_t two_pad_max=4;
 441         Short_t two_row_max=3;
 442
 443         if(abs(max_row - (*mp2)->GetRow())<two_row_max &&
 444            abs(max_np - (*mp2)->GetPad())<two_pad_max  &&
 445            abs(max_nt - (*mp2)->GetTime())<7){
 446
 447           cluster_value+=(*mp2)->GetSum();
 448           n_used_pads++;
 449
 450           AliSimDigits *digrow_tmp1;
 451           if(fRawData){
 452             digrow_tmp1 = (AliSimDigits*)fDigarr->GetRow(sec,(*mp2)->GetRow());
 453           }else{
 454             digrow_tmp1 = (AliSimDigits*)fDigarr->LoadRow(sec,(*mp2)->GetRow());
 455           }
 456           digrow_tmp1->ExpandBuffer(); //decrunch
 457
 458           for(Short_t tb=(*mp2)->GetBegin(); tb<(*mp2)->GetEnd()+1; tb++){
 459             Short_t adc_tmp = digrow_tmp1->GetDigitFast(tb,(*mp2)->GetPad());
 460             AliTPCvtpr *vtpr=new AliTPCvtpr(adc_tmp,tb,(*mp2)->GetPad(),(*mp2)->GetRow());
 461               tmp->fCluster.push_back(vtpr);
 462           }
 463
 464           //which one is bigger
 465           if( (*mp2)->GetAdc() > max_dig ){
 466             max_dig      =(*mp2)->GetAdc() ;
 467             max_sum_adc  =(*mp2)->GetSum() ;
 468             max_nt       =(*mp2)->GetTime();
 469             max_np       =(*mp2)->GetPad() ;
 470             max_row      =(*mp2)->GetRow() ;
 471           } else if ( (*mp2)->GetAdc() == max_dig ){
 472             if( (*mp2)->GetSum() > max_sum_adc){
 473               max_dig      =(*mp2)->GetAdc() ;
 474               max_sum_adc  =(*mp2)->GetSum() ;
 475               max_nt       =(*mp2)->GetTime();
 476               max_np       =(*mp2)->GetPad() ;
 477               max_row      =(*mp2)->GetRow() ;
 478             }
 479           }
 480           maxima_in_sector.erase(mp2);
 481           mp2--;
 482         }
 483       }//inside loop
 484
 485       //through out ADC=6,7 on 1 pad, 2 tb and ADC=12 on 2 pads,2 tb
 486       if(n_used_pads==1 && cluster_value/tmp->fCluster.size()<3.6)continue;
 487       if(n_used_pads==2 && cluster_value/tmp->fCluster.size()<3.1)continue;
 488
 489       tmp->SetADCcluster(cluster_value);
 490       tmp->SetNpads(n_used_pads);
 491       tmp->SetMax(AliTPCvtpr(max_dig,max_nt,max_np,max_row));
 492       tmp->SetSec(sec);
 493       tmp->SetSize();
 494
 495       cluster_counter++;
 496
 497       //save each cluster into file
 498
 499       AliTPCClustersRow *clrow= new AliTPCClustersRow();
 500       fRowCl=clrow;
 501       clrow->SetClass("AliTPCclusterKr");
 502       clrow->SetArray(1);
 503       fOutput->GetBranch("Segment")->SetAddress(&clrow);
 504
 505       Int_t tmp_cluster=0;
 506       TClonesArray * arr = fRowCl->GetArray();
 507       AliTPCclusterKr * cl = new ((*arr)[tmp_cluster]) AliTPCclusterKr(*tmp);
 508
 509       fOutput->Fill();
 510       delete clrow;
 511
 512       //end of save each cluster into file adc.root
 513     }//outer loop
 514
 515   }//end sector for
 516   cout<<"Number of clusters in event: "<<cluster_counter<<endl;
 517
 518   return 0;
 519 }
 520
 521 ////____________________________________________________________________________