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