TPC/Sim/AliTPCDigitizer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2000, 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$ */
  17
  18 /*
  19   Class for creating of the sumable digits and digits from MC data
  20   //
  21   The input :  ideal signals (Hits->Diffusion->Attachment -Ideal signal)
  22   The output:  raw digits
  23
  24   Effect implemented:
  25   1. Pad by pad gain map
  26   2. Noise map
  27   3. The dead channels identified  - zerro noise for corresponding pads
  28      In this case the outpu equal zerro
  29
  30 */
  31
  32
  33
  34
  35 #include <stdlib.h>
  36 #include <TTree.h>
  37 #include <TObjArray.h>
  38 #include <TFile.h>
  39 #include <TDirectory.h>
  40 #include <Riostream.h>
  41 #include <TParameter.h>
  42
  43 #include "AliTPCDigitizer.h"
  44
  45 #include "AliTPC.h"
  46 #include "AliTPCParam.h"
  47 #include "AliTPCParamSR.h"
  48 #include "AliRun.h"
  49 #include "AliLoader.h"
  50 #include "AliPDG.h"
  51 #include "AliDigitizationInput.h"
  52 #include "AliSimDigits.h"
  53 #include "AliLog.h"
  54
  55 #include "AliTPCcalibDB.h"
  56 #include "AliTPCCalPad.h"
  57 #include "AliTPCCalROC.h"
  58
  59 using std::cout;
  60 using std::cerr;
  61 using std::endl;
  62 ClassImp(AliTPCDigitizer)
  63
  64 //___________________________________________
  65   AliTPCDigitizer::AliTPCDigitizer() :AliDigitizer(),fDebug(0)
  66 {
  67   //
  68 // Default ctor - don't use it
  69 //
  70
  71 }
  72
  73 //___________________________________________
  74 AliTPCDigitizer::AliTPCDigitizer(AliDigitizationInput* digInput)
  75   :AliDigitizer(digInput),fDebug(0)
  76 {
  77   //
  78 // ctor which should be used
  79 //
  80   AliDebug(2,"(AliDigitizationInput* digInput) was processed");
  81 }
  82
  83 //------------------------------------------------------------------------
  84 AliTPCDigitizer::~AliTPCDigitizer()
  85 {
  86 // Destructor
  87 }
  88
  89
  90
  91 //------------------------------------------------------------------------
  92 Bool_t AliTPCDigitizer::Init()
  93 {
  94 // Initialization
  95
  96  return kTRUE;
  97 }
  98
  99
 100 //------------------------------------------------------------------------
 101 void AliTPCDigitizer::Digitize(Option_t* option)
 102 {
 103   DigitizeFast(option);
 104 }
 105 //------------------------------------------------------------------------
 106 void AliTPCDigitizer::DigitizeFast(Option_t* option)
 107 {
 108
 109   // merge input tree's with summable digits
 110   //output stored in TreeTPCD
 111   char s[100];
 112   char ss[100];
 113   TString optionString = option;
 114   if (!strcmp(optionString.Data(),"deb")) {
 115     cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
 116     fDebug = 3;
 117   }
 118   //get detector and geometry
 119
 120
 121   AliRunLoader *rl, *orl;
 122   AliLoader *gime, *ogime;
 123
 124   if (gAlice == 0x0)
 125    {
 126      Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
 127      rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 128      if (rl == 0x0)
 129       {
 130         Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
 131         return;
 132       }
 133      rl->LoadgAlice();
 134      rl->GetAliRun();
 135    }
 136   AliTPC *pTPC  = (AliTPC *) gAlice->GetModule("TPC");
 137   AliTPCParam * param = pTPC->GetParam();
 138
 139   //sprintf(s,param->GetTitle());
 140   snprintf(s,100,"%s",param->GetTitle());
 141   //sprintf(ss,"75x40_100x60");
 142   snprintf(ss,100,"75x40_100x60");
 143   if(strcmp(s,ss)==0){
 144     printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
 145     delete param;
 146     param=new AliTPCParamSR();
 147   }
 148   else{
 149     //sprintf(ss,"75x40_100x60_150x60");
 150     snprintf(ss,100,"75x40_100x60_150x60");
 151    if(strcmp(s,ss)!=0) {
 152      printf("No TPC parameters found...\n");
 153      exit(2);
 154    }
 155   }
 156
 157   pTPC->GenerNoise(500000); //create table with noise
 158   //
 159   Int_t nInputs = fDigInput->GetNinputs();
 160   Int_t * masks = new Int_t[nInputs];
 161   for (Int_t i=0; i<nInputs;i++)
 162     masks[i]= fDigInput->GetMask(i);
 163   Short_t **pdig= new Short_t*[nInputs];   //pointers to the expanded digits array
 164   Int_t **ptr=  new Int_t*[nInputs];       //pointers to the expanded tracks array
 165   Bool_t *active=  new Bool_t[nInputs];    //flag for active input segments
 166   Char_t phname[100];
 167
 168   //create digits array for given sectors
 169   // make indexes
 170   //
 171   //create branch's in TPC treeD
 172   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 173   ogime = orl->GetLoader("TPCLoader");
 174   TTree * tree  = ogime->TreeD();
 175   AliSimDigits * digrow = new AliSimDigits;
 176
 177   if (tree == 0x0)
 178    {
 179      ogime->MakeTree("D");
 180      tree  = ogime->TreeD();
 181    }
 182   tree->Branch("Segment","AliSimDigits",&digrow);
 183   //
 184   AliSimDigits ** digarr = new AliSimDigits*[nInputs];
 185   for (Int_t i1=0;i1<nInputs; i1++)
 186     {
 187       digarr[i1]=0;
 188      //    intree[i1]
 189       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
 190       gime = rl->GetLoader("TPCLoader");
 191       gime->LoadSDigits("read");
 192       TTree * treear =  gime->TreeS();
 193
 194       if (!treear)
 195        {
 196         cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
 197             <<" input "<< i1<<endl;
 198         for (Int_t i2=0;i2<i1+1; i2++){
 199
 200           if(digarr[i2])  delete digarr[i2];
 201         }
 202         delete [] digarr;
 203         delete [] active;
 204         delete []masks;
 205         delete []pdig;
 206         delete []ptr;
 207         return;
 208        }
 209
 210       //sprintf(phname,"lhcphase%d",i1);
 211       snprintf(phname,100,"lhcphase%d",i1);
 212       TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
 213                                ->FindObject("lhcphase0");
 214       if(!ph){
 215         cerr<<"AliTPCDigitizer: LHC phase  not found in"
 216             <<" input "<< i1<<endl;
 217         for (Int_t i2=0;i2<i1+1; i2++){
 218           if(digarr[i2])  delete digarr[i2];
 219         }
 220         delete [] digarr;
 221         delete [] active;
 222         delete []masks;
 223         delete []pdig;
 224         delete []ptr;
 225         return;
 226       }
 227       tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
 228               //
 229       if (treear->GetIndex()==0)
 230         treear->BuildIndex("fSegmentID","fSegmentID");
 231       treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
 232     }
 233
 234
 235
 236
 237   //
 238
 239   param->SetZeroSup(2);
 240
 241   Int_t zerosup = param->GetZeroSup();
 242   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 243   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 244   //
 245   //Loop over segments of the TPC
 246
 247   for (Int_t segmentID=0; segmentID<param->GetNRowsTotal(); segmentID++)
 248    {
 249     Int_t sec, row;
 250     if (!param->AdjustSectorRow(segmentID,sec,row))
 251      {
 252       cerr<<"AliTPC warning: invalid segment ID ! "<<segmentID<<endl;
 253       continue;
 254      }
 255     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sec);  // pad gains per given sector
 256     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sec);  // noise per given sector
 257     digrow->SetID(segmentID);
 258
 259     Int_t nrows = 0;
 260     Int_t ncols = 0;
 261
 262     Bool_t digitize = kFALSE;
 263     for (Int_t i=0;i<nInputs; i++)
 264      {
 265
 266       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 267       gime = rl->GetLoader("TPCLoader");
 268
 269       if (gime->TreeS()->GetEntryWithIndex(segmentID,segmentID) >= 0) {
 270         digarr[i]->ExpandBuffer();
 271         digarr[i]->ExpandTrackBuffer();
 272         nrows = digarr[i]->GetNRows();
 273         ncols = digarr[i]->GetNCols();
 274         active[i] = kTRUE;
 275         if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
 276       } else {
 277         active[i] = kFALSE;
 278       }
 279       if (GetRegionOfInterest() && !digitize) break;
 280      }
 281     if (!digitize) continue;
 282
 283     digrow->Allocate(nrows,ncols);
 284     digrow->AllocateTrack(3);
 285
 286     Float_t q=0;
 287     Int_t label[1000]; //stack for 300 events
 288     Int_t labptr = 0;
 289
 290     Int_t nElems = nrows*ncols;
 291
 292     for (Int_t i=0;i<nInputs; i++)
 293      if (active[i]) {
 294        pdig[i] = digarr[i]->GetDigits();
 295        ptr[i]  = digarr[i]->GetTracks();
 296       }
 297
 298     Short_t *pdig1= digrow->GetDigits();
 299     Int_t   *ptr1= digrow->GetTracks() ;
 300
 301
 302
 303     for (Int_t elem=0;elem<nElems; elem++)
 304      {
 305
 306        q=0;
 307        labptr=0;
 308        // looop over digits
 309         for (Int_t i=0;i<nInputs; i++) if (active[i])
 310          {
 311           //          q  += digarr[i]->GetDigitFast(rows,col);
 312             q  += *(pdig[i]);
 313
 314            for (Int_t tr=0;tr<3;tr++)
 315             {
 316              //             Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
 317              Int_t lab = ptr[i][tr*nElems];
 318              if ( (lab > 1) && *(pdig[i])>zerosup)
 319               {
 320                 label[labptr]=lab+masks[i];
 321                 labptr++;
 322               }
 323             }
 324            pdig[i]++;
 325            ptr[i]++;
 326          }
 327         q/=16.;  //conversion factor
 328         Float_t gain = gainROC->GetValue(row,elem/nrows);  // get gain for given - pad-row pad
 329         //if (gain<0.5){
 330           //printf("problem\n");
 331         //}
 332         q*= gain;
 333         Float_t noisePad = noiseROC->GetValue(row,elem/nrows);
 334         //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
 335         Float_t noise  = pTPC->GetNoise();
 336         q+=noise*noisePad;
 337         q=TMath::Nint(q);
 338         if (q > zerosup)
 339          {
 340           if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
 341           //digrow->SetDigitFast((Short_t)q,rows,col);
 342           *pdig1 =Short_t(q);
 343           for (Int_t tr=0;tr<3;tr++)
 344            {
 345             if (tr<labptr)
 346              ptr1[tr*nElems] = label[tr];
 347            }
 348           }
 349         pdig1++;
 350         ptr1++;
 351      }
 352     //
 353     //  glitch filter
 354     //
 355     digrow->GlitchFilter();
 356     //
 357     digrow->CompresBuffer(1,zerosup);
 358     digrow->CompresTrackBuffer(1);
 359     tree->Fill();
 360     if (fDebug>0) cerr<<sec<<"\t"<<row<<"\n";
 361    } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
 362
 363
 364   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 365   ogime = orl->GetLoader("TPCLoader");
 366   ogime->WriteDigits("OVERWRITE");
 367
 368   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
 369
 370   delete digrow;
 371   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
 372   delete []masks;
 373   delete []pdig;
 374   delete []ptr;
 375   delete []active;
 376   delete []digarr;
 377 }
 378
 379
 380
 381 //------------------------------------------------------------------------
 382 void AliTPCDigitizer::DigitizeSave(Option_t* option)
 383 {
 384
 385   // merge input tree's with summable digits
 386   //output stored in TreeTPCD
 387
 388   TString optionString = option;
 389   if (!strcmp(optionString.Data(),"deb")) {
 390     cout<<"AliTPCDigitizer::Digitize: called with option deb "<<endl;
 391     fDebug = 3;
 392   }
 393   //get detector and geometry
 394   AliRunLoader *rl, *orl;
 395   AliLoader *gime, *ogime;
 396
 397
 398   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 399   ogime = orl->GetLoader("TPCLoader");
 400
 401   rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 402   //gime = rl->GetLoader("TPCLoader");
 403   rl->GetLoader("TPCLoader");
 404   rl->LoadgAlice();
 405   AliRun* alirun = rl->GetAliRun();
 406
 407   AliTPC *pTPC  = (AliTPC *) alirun->GetModule("TPC");
 408   AliTPCParam * param = pTPC->GetParam();
 409   pTPC->GenerNoise(500000); //create teble with noise
 410   printf("noise %f \n",  param->GetNoise()*param->GetNoiseNormFac());
 411   //
 412   Int_t nInputs = fDigInput->GetNinputs();
 413   // stupid protection...
 414   if (nInputs <= 0) return;
 415   //
 416   Int_t * masks = new Int_t[nInputs];
 417   for (Int_t i=0; i<nInputs;i++)
 418     masks[i]= fDigInput->GetMask(i);
 419
 420   AliSimDigits ** digarr = new AliSimDigits*[nInputs];
 421   for(Int_t ii=0;ii<nInputs;ii++) digarr[ii]=0;
 422
 423   for (Int_t i1=0;i1<nInputs; i1++)
 424    {
 425      //digarr[i1]=0;
 426     //    intree[i1]
 427     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
 428     gime = rl->GetLoader("TPCLoader");
 429
 430     TTree * treear =  gime->TreeS();
 431     //
 432     if (!treear) {
 433       cerr<<" TPC -  not existing input = \n"<<i1<<" ";
 434       delete [] masks;
 435       for(Int_t i=0; i<nInputs; i++) delete digarr[i];
 436       delete [] digarr;
 437       return;
 438     }
 439     //
 440     TBranch * br = treear->GetBranch("fSegmentID");
 441     if (br) br->GetFile()->cd();
 442     treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
 443    }
 444
 445   rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 446   gime = rl->GetLoader("TPCLoader");
 447   Stat_t nentries = gime->TreeS()->GetEntries();
 448
 449
 450   //create branch's in TPC treeD
 451   AliSimDigits * digrow = new AliSimDigits;
 452   TTree * tree  = ogime->TreeD();
 453
 454   tree->Branch("Segment","AliSimDigits",&digrow);
 455   param->SetZeroSup(2);
 456
 457   Int_t zerosup = param->GetZeroSup();
 458   //Loop over segments of the TPC
 459
 460   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 461   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 462   for (Int_t n=0; n<nentries; n++) {
 463     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 464     gime = rl->GetLoader("TPCLoader");
 465     gime->TreeS()->GetEvent(n);
 466
 467     digarr[0]->ExpandBuffer();
 468     digarr[0]->ExpandTrackBuffer();
 469
 470
 471     for (Int_t i=1;i<nInputs; i++){
 472 //      fDigInput->GetInputTreeTPCS(i)->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
 473       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 474       gime = rl->GetLoader("TPCLoader");
 475       gime->TreeS()->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
 476       digarr[i]->ExpandBuffer();
 477       digarr[i]->ExpandTrackBuffer();
 478       if ((digarr[0]->GetID()-digarr[i]->GetID())>0)
 479        printf("problem\n");
 480
 481     }
 482
 483     Int_t sec, row;
 484     if (!param->AdjustSectorRow(digarr[0]->GetID(),sec,row)) {
 485       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
 486       continue;
 487     }
 488
 489     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sec);  // pad gains per given sector
 490     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sec);  // noise per given sector
 491     digrow->SetID(digarr[0]->GetID());
 492
 493     Int_t nrows = digarr[0]->GetNRows();
 494     Int_t ncols = digarr[0]->GetNCols();
 495     digrow->Allocate(nrows,ncols);
 496     digrow->AllocateTrack(3);
 497
 498     Float_t q=0;
 499     Int_t label[1000]; //stack for 300 events
 500     Int_t labptr = 0;
 501
 502
 503
 504     for (Int_t rows=0;rows<nrows; rows++){
 505       for (Int_t col=0;col<ncols; col++){
 506
 507        q=0;
 508        labptr=0;
 509        // looop over digits
 510         for (Int_t i=0;i<nInputs; i++){
 511          q  += digarr[i]->GetDigitFast(rows,col);
 512           //q  += *(pdig[i]);
 513
 514           for (Int_t tr=0;tr<3;tr++) {
 515            Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
 516            //Int_t lab = ptr[i][tr*nElems];
 517             if ( (lab > 1) ) {
 518               label[labptr]=lab+masks[i];
 519               labptr++;
 520             }
 521           }
 522          // pdig[i]++;
 523          //ptr[i]++;
 524
 525         }
 526        q/=16.;  //conversion factor
 527        //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
 528        Float_t gain = gainROC->GetValue(row,col);
 529        q*= gain;
 530        Float_t noisePad = noiseROC->GetValue(row, col);
 531
 532        Float_t noise  = pTPC->GetNoise();
 533        q+=noise*noisePad;
 534
 535         q=TMath::Nint(q);
 536         if (q > zerosup){
 537
 538          if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
 539          digrow->SetDigitFast((Short_t)q,rows,col);
 540          // *pdig1 =Short_t(q);
 541          for (Int_t tr=0;tr<3;tr++){
 542            if (tr<labptr)
 543              ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],rows,col,tr);
 544            //ptr1[tr*nElems] = label[tr];
 545            //else
 546              //           ((AliSimDigits*)digrow)->SetTrackIDFast(-1,rows,col,tr);
 547            //  ptr1[tr*nElems] = 1;
 548          }
 549        }
 550        //pdig1++;
 551        //ptr1++;
 552     }
 553     }
 554
 555     digrow->CompresBuffer(1,zerosup);
 556     digrow->CompresTrackBuffer(1);
 557     tree->Fill();
 558     if (fDebug>0) cerr<<sec<<"\t"<<row<<"\n";
 559   }
 560 //  printf("end TPC merging - end -Tree %s\t%p\n",fDigInput->GetInputTreeH(0)->GetName(),fDigInput->GetInputTreeH(0)->GetListOfBranches()->At(3));
 561   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
 562     ogime->WriteDigits("OVERWRITE");
 563
 564     for (Int_t i=1;i<nInputs; i++)
 565      {
 566       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 567       gime = rl->GetLoader("TPCLoader");
 568       gime->UnloadSDigits();
 569      }
 570     ogime->UnloadDigits();
 571
 572   delete digrow;
 573   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
 574   delete [] masks;
 575   delete [] digarr;
 576 }