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