T0/AliT0Reconstructor.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$ */
  17 /*********************************************************************
  18  *  T0 reconstruction and filling ESD
  19  *  - reconstruct mean time (interation time)
  20  *  - vertex position
  21  *  -  multiplicity
  22  ********************************************************************/
  23
  24 #include <AliESDEvent.h>
  25 #include "AliLog.h"
  26 #include "AliT0RecPoint.h"
  27 #include "AliRawReader.h"
  28 #include "AliT0RawReader.h"
  29 #include "AliT0digit.h"
  30 #include "AliT0Reconstructor.h"
  31 #include "AliT0Parameters.h"
  32 #include "AliT0Calibrator.h"
  33
  34 #include <TArrayI.h>
  35 #include <TGraph.h>
  36 #include <TMath.h>
  37 #include <iostream.h>
  38
  39 ClassImp(AliT0Reconstructor)
  40
  41   AliT0Reconstructor:: AliT0Reconstructor(): AliReconstructor(),
  42                                              fdZonA(0),
  43                                              fdZonC(0),
  44                                              fZposition(0),
  45                                              fParam(NULL),
  46                                              fAmpLEDrec(),
  47                                              fCalib()
  48 {
  49   //constructor
  50
  51  AliDebug(1,"Start reconstructor ");
  52
  53   fParam = AliT0Parameters::Instance();
  54   fParam->Init();
  55
  56   for (Int_t i=0; i<24; i++){
  57         TGraph* gr = fParam ->GetAmpLEDRec(i);
  58         if (gr) fAmpLEDrec.AddAtAndExpand(gr,i) ;
  59 }
  60
  61   fdZonC = TMath::Abs(fParam->GetZPositionShift("T0/C/PMT1"));
  62   fdZonA = TMath::Abs(fParam->GetZPositionShift("T0/A/PMT15"));
  63   fCalib = new AliT0Calibrator();
  64 }
  65 //____________________________________________________________________
  66
  67 AliT0Reconstructor::AliT0Reconstructor(const AliT0Reconstructor &r):
  68   AliReconstructor(r),
  69                                              fdZonA(0),
  70                                              fdZonC(0),
  71                                              fZposition(0),
  72                                              fParam(NULL),
  73                                              fAmpLEDrec(),
  74                                              fCalib()
  75
  76  {
  77   //
  78   // AliT0Reconstructor copy constructor
  79   //
  80
  81   ((AliT0Reconstructor &) r).Copy(*this);
  82
  83 }
  84
  85 //_____________________________________________________________________________
  86 AliT0Reconstructor &AliT0Reconstructor::operator=(const AliT0Reconstructor &r)
  87 {
  88   //
  89   // Assignment operator
  90   //
  91
  92   if (this != &r) ((AliT0Reconstructor &) r).Copy(*this);
  93   return *this;
  94
  95 }
  96
  97 //_____________________________________________________________________________
  98
  99 void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const
 100
 101 {
 102   // T0 digits reconstruction
 103   // T0RecPoint writing
 104
 105
 106   TArrayI * timeCFD = new TArrayI(24);
 107   TArrayI * timeLED = new TArrayI(24);
 108   TArrayI * chargeQT0 = new TArrayI(24);
 109   TArrayI * chargeQT1 = new TArrayI(24);
 110
 111
 112   //  Int_t mV2Mip = param->GetmV2Mip();
 113   //mV2Mip = param->GetmV2Mip();
 114   Float_t channelWidth = fParam->GetChannelWidth() ;
 115   Int_t meanT0 = fParam->GetMeanT0();
 116
 117   AliDebug(1,Form("Start DIGITS reconstruction "));
 118
 119
 120  TBranch *brDigits=digitsTree->GetBranch("T0");
 121   AliT0digit *fDigits = new AliT0digit() ;
 122   if (brDigits) {
 123     brDigits->SetAddress(&fDigits);
 124   }else{
 125     AliError(Form("EXEC Branch T0 digits not found"));
 126      return;
 127   }
 128
 129   digitsTree->GetEvent(0);
 130   digitsTree->GetEntry(0);
 131   brDigits->GetEntry(0);
 132   fDigits->GetTimeCFD(*timeCFD);
 133   fDigits->GetTimeLED(*timeLED);
 134   fDigits->GetQT0(*chargeQT0);
 135   fDigits->GetQT1(*chargeQT1);
 136
 137
 138   Float_t besttimeA=999999;
 139   Float_t besttimeC=999999;
 140   Int_t pmtBestA=99999;
 141   Int_t pmtBestC=99999;
 142   Float_t timeDiff=999999, meanTime=0;
 143
 144
 145   AliT0RecPoint* frecpoints= new AliT0RecPoint ();
 146   clustersTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
 147
 148   Float_t time[24], adc[24];
 149   for (Int_t ipmt=0; ipmt<24; ipmt++) {
 150     if(timeCFD->At(ipmt)>0 ){
 151       Double_t qt0 = Double_t(chargeQT0->At(ipmt));
 152       Double_t qt1 = Double_t(chargeQT1->At(ipmt));
 153       if((qt1-qt0)>0)  adc[ipmt] = TMath::Exp( Double_t (channelWidth*(qt1-qt0)/1000));
 154       time[ipmt] = fCalib-> WalkCorrection( ipmt, Int_t(qt1) , timeCFD->At(ipmt), "pdc" ) ;
 155
 156       //LED
 157       Double_t sl = (timeLED->At(ipmt) - time[ipmt])*channelWidth;
 158       Double_t qt=((TGraph*)fAmpLEDrec.At(ipmt))->Eval(sl/1000.);
 159       frecpoints->SetTime(ipmt,time[ipmt]);
 160       frecpoints->SetAmp(ipmt,adc[ipmt]);
 161       frecpoints->SetAmpLED(ipmt,qt);
 162     }
 163     else {
 164       time[ipmt] = 0;
 165       adc[ipmt] = 0;
 166     }
 167   }
 168
 169   for (Int_t ipmt=0; ipmt<12; ipmt++){
 170     if(time[ipmt] > 1 ) {
 171       if(time[ipmt]<besttimeC){
 172         besttimeC=time[ipmt]; //timeC
 173         pmtBestC=ipmt;
 174       }
 175     }
 176   }
 177   for ( Int_t ipmt=12; ipmt<24; ipmt++){
 178     if(time[ipmt] > 1) {
 179       if(time[ipmt]<besttimeA) {
 180         besttimeA=time[ipmt]; //timeA
 181         pmtBestA=ipmt;}
 182     }
 183   }
 184   if(besttimeA !=999999)  frecpoints->SetTimeBestA(Int_t(besttimeA));
 185   if( besttimeC != 999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
 186   AliDebug(1,Form(" besttimeA %f ps,  besttimeC %f ps",besttimeA, besttimeC));
 187   Float_t c = 0.0299792; // cm/ps
 188   Float_t vertex = 0;
 189   if(besttimeA !=999999 && besttimeC != 999999 ){
 190     timeDiff =(besttimeC - besttimeA)*channelWidth;
 191     meanTime = (meanT0 - (besttimeA + besttimeC)/2) * channelWidth;
 192     vertex = c*(timeDiff)/2. + (fdZonA - fdZonC)/2; //-(lenr-lenl))/2;
 193     AliDebug(1,Form("  timeDiff %f ps,  meanTime %f ps, vertex %f cm",timeDiff, meanTime,vertex ));
 194     frecpoints->SetVertex(vertex);
 195     frecpoints->SetMeanTime(Int_t(meanTime));
 196
 197   }
 198   //time in each channel as time[ipmt]-MeanTimeinThisChannel(with vertex=0)
 199   for (Int_t ipmt=0; ipmt<24; ipmt++) {
 200     if(time[ipmt]>1) {
 201 //      time[ipmt] = (time[ipmt] - fTime0vertex[ipmt])*channelWidth;
 202       time[ipmt] = time[ipmt] * channelWidth;
 203       frecpoints->SetTime(ipmt,time[ipmt]);
 204     }
 205   }
 206   clustersTree->Fill();
 207
 208   delete timeCFD;
 209   delete timeLED;
 210   delete chargeQT0;
 211   delete chargeQT1;
 212 }
 213
 214
 215 //_______________________________________________________________________
 216
 217 void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
 218 {
 219   // T0 raw ->
 220   // T0RecPoint writing
 221
 222   //Q->T-> coefficients !!!! should be measured!!!
 223   Int_t allData[110][5];
 224
 225   Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
 226   TString option = GetOption();
 227    AliDebug(10,Form("Option: %s\n", option.Data()));
 228
 229    for (Int_t i0=0; i0<105; i0++)
 230      {
 231        for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0;
 232      }
 233
 234    Int_t besttimeA=9999999;
 235    Int_t besttimeC=9999999;
 236    Int_t pmtBestA=99999;
 237    Int_t pmtBestC=99999;
 238    Int_t timeDiff=9999999, meanTime=0;
 239    Double_t qt=0;
 240
 241    AliT0RecPoint* frecpoints= new AliT0RecPoint ();
 242
 243    recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
 244
 245
 246    AliDebug(10," before read data ");
 247    AliT0RawReader myrawreader(rawReader);
 248    if (!myrawreader.Next())
 249      AliDebug(1,Form(" no raw data found!!"));
 250    else
 251      {
 252        for (Int_t i=0; i<105; i++) {
 253          for (Int_t iHit=0; iHit<5; iHit++)
 254            {
 255              allData[i][iHit] = myrawreader.GetData(i,iHit);
 256            }
 257        }
 258
 259        Float_t channelWidth = fParam->GetChannelWidth() ;
 260
 261        Int_t meanT0 = fParam->GetMeanT0();
 262        if(option == "pdc"){
 263          for (Int_t in=0; in<24; in++)
 264            {
 265
 266              timeLED[in] = allData[in+1][0] ;
 267              timeCFD[in] = allData[in+25][0] ;
 268              chargeQT1[in] = allData[in+57][0] ;
 269              chargeQT0[in] = allData[in+80][0] ;
 270            }
 271        }
 272
 273        if(option == "cosmic") {
 274          for (Int_t in=0; in<12; in++)
 275            {
 276              timeCFD[in] = allData[in+1][0] ;
 277              timeCFD[in+12] = allData[in+56+1][0] ;
 278              timeLED[in] = allData[in+12+1][0] ;
 279              timeLED[in+12] = allData[in+68+1][0] ;
 280            }
 281
 282          for (Int_t in=0; in<24;  in=in+2)
 283            {
 284              Int_t cc=in/2;
 285              chargeQT1[cc]=allData[in+25][0];
 286              chargeQT0[cc]=allData[in+26][0];
 287            }
 288          for (Int_t in=24; in<48;  in=in+2)
 289            {
 290              Int_t cc=in/2;
 291              chargeQT1[cc]=allData[in+57][0];
 292              chargeQT0[cc]=allData[in+58][0];
 293            }
 294
 295        }
 296        for (Int_t in=0; in<24; in++)
 297          AliDebug(10, Form(" readed Raw %i %i %i %i %i", in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in]));
 298
 299
 300        Int_t time[24], adc[24];
 301        for (Int_t ipmt=0; ipmt<24; ipmt++) {
 302          if(timeCFD[ipmt]>0 && timeLED[ipmt]>0){
 303
 304            if(option == "pdc"){
 305              Double_t qt0 = Double_t(chargeQT0[ipmt]);
 306              Double_t qt1 = Double_t(chargeQT1[ipmt]);
 307              if((qt1-qt0)>0)  adc[ipmt] = Int_t(TMath::Exp( Double_t (channelWidth*(qt1-qt0)/1000.)));
 308              time[ipmt] = fCalib-> WalkCorrection( ipmt,Int_t(qt1) , timeCFD[ipmt], "pdc" ) ;
 309              Double_t sl = (timeLED[ipmt] - time[ipmt])*channelWidth;
 310              if(fAmpLEDrec.At(ipmt))
 311                qt=((TGraph*)fAmpLEDrec.At(ipmt))->Eval(sl/1000.);
 312              frecpoints->SetTime(ipmt,time[ipmt]);
 313              frecpoints->SetAmp(ipmt,adc[ipmt]);
 314              frecpoints->SetAmpLED(ipmt,qt);
 315              AliDebug(10,Form(" QTC %i mv,  time in chann %i ",adc[ipmt] ,time[ipmt]));
 316            }
 317            if(option == "cosmic") {
 318              //      if(ipmt == 15) continue; //skip crashed PMT
 319              if(( chargeQT1[ipmt] - chargeQT0[ipmt])>0)
 320                adc[ipmt] = chargeQT1[ipmt] - chargeQT0[ipmt];
 321              else
 322                adc[ipmt] = 0;
 323              //      time[ipmt] = fCalib-> WalkCorrection( ipmt, adc[ipmt], timeCFD[ipmt],"cosmic" ) ;
 324              // time[ipmt] =  timeCFD[ipmt] ;
 325              Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
 326              time[ipmt] = fCalib-> WalkCorrection( ipmt, Int_t(sl), timeCFD[ipmt],"cosmic" ) ;
 327              //      if(fAmpLEDrec.At(ipmt))
 328              // qt=((TGraph*)fAmpLEDrec.At(ipmt))->Eval(sl);
 329              time[ipmt] = time[ipmt] - allData[0][0] + 5000;
 330              AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
 331                               ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
 332              frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
 333              frecpoints->SetAmp(ipmt, Float_t(adc[ipmt]));
 334              frecpoints->SetAmpLED(ipmt, Float_t(qt));
 335            }
 336
 337          }
 338          else {
 339            time[ipmt] = 0;
 340            adc[ipmt] = 0;
 341          }
 342        }
 343
 344        for (Int_t ipmt=0; ipmt<12; ipmt++){
 345          if(time[ipmt] > 1 ) {
 346            if(time[ipmt]<besttimeC){
 347              besttimeC=time[ipmt]; //timeC
 348              pmtBestC=ipmt;
 349            }
 350          }
 351        }
 352        for ( Int_t ipmt=12; ipmt<24; ipmt++){
 353          if(time[ipmt] > 1) {
 354            if(time[ipmt]<besttimeA) {
 355              besttimeA=time[ipmt]; //timeA
 356              pmtBestA=ipmt;}
 357          }
 358        }
 359        if(besttimeA !=9999999)  frecpoints->SetTimeBestA(Int_t(besttimeA));
 360        if( besttimeC != 9999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
 361        AliDebug(1,Form(" besttimeA %i ps,  besttimeC %i ps",besttimeA, besttimeC));
 362        Float_t c = 0.0299792; // cm/ps
 363        Float_t vertex = 99999;
 364        if(besttimeA <9999999 && besttimeC < 9999999 ){
 365          timeDiff =Int_t (( besttimeC - besttimeA) *channelWidth);
 366          if(option == "pdc")
 367            meanTime = Int_t((meanT0 - (besttimeA + besttimeC)/2) * channelWidth);
 368          if(option == "cosmic")   meanTime =  (besttimeA + besttimeC)/2;
 369          vertex = c*(Float_t(timeDiff))/2.+ (fdZonA - fdZonC)/2;
 370          AliDebug(1,Form("  timeDiff %i ps,  meanTime %i ps, vertex %f cm",timeDiff, meanTime,vertex ));
 371          frecpoints->SetVertex(vertex);
 372          frecpoints->SetMeanTime(Int_t(meanTime));
 373
 374        }
 375      } // if (else )raw data
 376    recTree->Fill();
 377    if(frecpoints) delete frecpoints;
 378 }
 379 //____________________________________________________________
 380
 381 void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const
 382 {
 383
 384   /***************************************************
 385   Resonstruct digits to vertex position
 386   ****************************************************/
 387
 388   AliDebug(1,Form("Start FillESD T0"));
 389
 390   TTree *treeR = clustersTree;
 391
 392    AliT0RecPoint* frecpoints= new AliT0RecPoint ();
 393     if (!frecpoints) {
 394     AliError("Reconstruct Fill ESD >> no recpoints found");
 395     return;
 396   }
 397
 398   AliDebug(1,Form("Start FillESD T0"));
 399   TBranch *brRec = treeR->GetBranch("T0");
 400   if (brRec) {
 401     brRec->SetAddress(&frecpoints);
 402   }else{
 403     AliError(Form("EXEC Branch T0 rec not found"));
 404     return;
 405   }
 406
 407     brRec->GetEntry(0);
 408     Float_t amp[24], time[24];
 409     Float_t  zPosition = frecpoints -> GetVertex();
 410     Float_t timeStart = frecpoints -> GetMeanTime() ;
 411     for ( Int_t i=0; i<24; i++) {
 412       time[i] = Float_t (frecpoints -> GetTime(i)); // ps to ns
 413       amp[i] = frecpoints -> GetAmp(i);
 414     }
 415     pESD->SetT0zVertex(zPosition); //vertex Z position
 416     pESD->SetT0(timeStart);        // interaction time
 417     pESD->SetT0time(time);         // best TOF on each PMT
 418     pESD->SetT0amplitude(amp);     // number of particles(MIPs) on each PMT
 419
 420     AliDebug(1,Form(" Z position %f cm,  T0  %f ps",zPosition , timeStart));
 421
 422 } // vertex in 3 sigma
 423
 424
 425
 426
 427
 428