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