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