X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=T0%2FAliT0Reconstructor.cxx;h=aadb1e74e13e7d1f755ad52ac2004d7fac9e695c;hb=e73860407b65d47729d5c56690eae63f99968b5c;hp=0d7d7c6d5e2c0f268131d91ac7d0c4ff6dc20576;hpb=69ad950f17a5543cf30f7739a9d455542261b4a9;p=u%2Fmrichter%2FAliRoot.git diff --git a/T0/AliT0Reconstructor.cxx b/T0/AliT0Reconstructor.cxx index 0d7d7c6d5e2..aadb1e74e13 100644 --- a/T0/AliT0Reconstructor.cxx +++ b/T0/AliT0Reconstructor.cxx @@ -34,6 +34,7 @@ #include #include #include +#include ClassImp(AliT0Reconstructor) @@ -42,76 +43,56 @@ ClassImp(AliT0Reconstructor) fdZonC(0), fZposition(0), fParam(NULL), - fAmpLEDrec() + fAmpLEDrec(), + fQTC(0), + fAmpLED(0), + fCalib() { //constructor - AliDebug(1,"Start reconstructor "); - fParam = AliT0Parameters::Instance(); fParam->Init(); + for (Int_t i=0; i<24; i++){ - TGraph* gr = fParam ->GetAmpLEDRec(i); - fAmpLEDrec.AddAtAndExpand(gr,i) ; -// fTime0vertex[i]= fParam->GetTimeV0(i); + TGraph* gr = fParam ->GetAmpLEDRec(i); + if (gr) fAmpLEDrec.AddAtAndExpand(gr,i) ; + TGraph* gr1 = fParam ->GetAmpLED(i); + if (gr1) fAmpLED.AddAtAndExpand(gr1,i) ; + TGraph* gr2 = fParam ->GetQTC(i); + if (gr2) fQTC.AddAtAndExpand(gr2,i) ; } + + fdZonC = TMath::Abs(fParam->GetZPositionShift("T0/C/PMT1")); fdZonA = TMath::Abs(fParam->GetZPositionShift("T0/A/PMT15")); -} -//____________________________________________________________________ - -AliT0Reconstructor::AliT0Reconstructor(const AliT0Reconstructor &r): - AliReconstructor(r), - fdZonA(0), - fdZonC(0), - fZposition(0), - fParam(NULL), - fAmpLEDrec() - { - // - // AliT0Reconstructor copy constructor - // - - ((AliT0Reconstructor &) r).Copy(*this); - -} - -//_____________________________________________________________________________ -AliT0Reconstructor &AliT0Reconstructor::operator=(const AliT0Reconstructor &r) -{ - // - // Assignment operator - // - - if (this != &r) ((AliT0Reconstructor &) r).Copy(*this); - return *this; + fCalib = new AliT0Calibrator(); } //_____________________________________________________________________________ - void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const { // T0 digits reconstruction - // T0RecPoint writing + Int_t refAmp = GetRecoParam()->GetRefAmp(); - TArrayI * timeCFD = new TArrayI(24); TArrayI * timeLED = new TArrayI(24); TArrayI * chargeQT0 = new TArrayI(24); TArrayI * chargeQT1 = new TArrayI(24); - AliT0Calibrator *calib=new AliT0Calibrator(); - - // Int_t mV2Mip = param->GetmV2Mip(); - //mV2Mip = param->GetmV2Mip(); + Float_t channelWidth = fParam->GetChannelWidth() ; - Int_t meanT0 = fParam->GetMeanT0(); + Float_t meanVertex = fParam->GetMeanVertex(); + Float_t c = 0.0299792; // cm/ps + Float_t vertex = 9999999; + Int_t timeDiff=999999, meanTime=999999, timeclock=999999; + AliDebug(1,Form("Start DIGITS reconstruction ")); + TBranch *brDigits=digitsTree->GetBranch("T0"); AliT0digit *fDigits = new AliT0digit() ; if (brDigits) { @@ -128,33 +109,41 @@ void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const fDigits->GetTimeLED(*timeLED); fDigits->GetQT0(*chargeQT0); fDigits->GetQT1(*chargeQT1); + Int_t onlineMean = fDigits->MeanTime(); + Int_t ref = fDigits->RefPoint(); Float_t besttimeA=999999; Float_t besttimeC=999999; Int_t pmtBestA=99999; Int_t pmtBestC=99999; - Float_t timeDiff=999999, meanTime=0; - - AliT0RecPoint* frecpoints= new AliT0RecPoint (); clustersTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1); Float_t time[24], adc[24]; for (Int_t ipmt=0; ipmt<24; ipmt++) { if(timeCFD->At(ipmt)>0 ){ - Double_t qt0 = Double_t(chargeQT0->At(ipmt)); - Double_t qt1 = Double_t(chargeQT1->At(ipmt)); - if((qt1-qt0)>0) adc[ipmt] = TMath::Exp( Double_t (channelWidth*(qt1-qt0)/1000)); - time[ipmt] = calib-> WalkCorrection( ipmt,Int_t(qt1) , timeCFD->At(ipmt), "pdc" ) ; + if(( chargeQT1->At(ipmt) - chargeQT0->At(ipmt))>0) + adc[ipmt] = chargeQT1->At(ipmt) - chargeQT0->At(ipmt); + else + adc[ipmt] = 0; + + time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD->At(ipmt)) ; + + Double_t sl = Double_t(timeLED->At(ipmt) - timeCFD->At(ipmt)); + // time[ipmt] = fCalib-> WalkCorrection( ipmt, Int_t(sl), timeCFD[ipmt],"cosmic" ) ; + AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ", + ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl))); + + Double_t ampMip =((TGraph*)fAmpLED.At(ipmt))->Eval(sl); + Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]); + AliDebug(10,Form(" Amlitude in MIPS LED %f , QTC %f in channels %i\n ",ampMip,qtMip, adc[ipmt])); + + frecpoints->SetTime(ipmt, Float_t(time[ipmt]) ); + frecpoints->SetAmp(ipmt, Float_t( ampMip)); //for cosmic &pp beam + frecpoints->SetAmpLED(ipmt, Float_t(qtMip)); - //LED - Double_t sl = (timeLED->At(ipmt) - time[ipmt])*channelWidth; - Double_t qt=((TGraph*)fAmpLEDrec.At(ipmt))->Eval(sl/1000.); - frecpoints->SetTime(ipmt,time[ipmt]); - frecpoints->SetAmp(ipmt,adc[ipmt]); - frecpoints->SetAmpLED(ipmt,qt); } else { time[ipmt] = 0; @@ -179,26 +168,23 @@ void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const } if(besttimeA !=999999) frecpoints->SetTimeBestA(Int_t(besttimeA)); if( besttimeC != 999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC)); - AliDebug(1,Form(" besttimeA %f ps, besttimeC %f ps",besttimeA, besttimeC)); - Float_t c = 0.0299792; // cm/ps - Float_t vertex = 0; - if(besttimeA !=999999 && besttimeC != 999999 ){ - timeDiff =(besttimeC - besttimeA)*channelWidth; - meanTime = (meanT0 - (besttimeA + besttimeC)/2) * channelWidth; - vertex = c*(timeDiff)/2. + (fdZonA - fdZonC)/2; //-(lenr-lenl))/2; - AliDebug(1,Form(" timeDiff %f ps, meanTime %f ps, vertex %f cm",timeDiff, meanTime,vertex )); + AliDebug(10,Form(" besttimeA %f ch, besttimeC %f ch",besttimeA, besttimeC)); + if(besttimeA !=999999 && besttimeC != 999999 ){ + // timeDiff = (besttimeC - besttimeA)*channelWidth; + timeDiff = (besttimeA - besttimeC)*channelWidth; + meanTime = Float_t((besttimeA + besttimeC)/2);// * channelWidth); + timeclock = Float_t(meanTime - ref); + vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2; + } frecpoints->SetVertex(vertex); - frecpoints->SetMeanTime(Int_t(meanTime)); + frecpoints->SetMeanTime(meanTime); + frecpoints->SetT0clock(timeclock); + //online mean + frecpoints->SetOnlineMean(Int_t(onlineMean)); + AliDebug(10,Form(" timeDiff %i #channel, meanTime %i #channel, vertex %f cm online mean %i timeclock %i ps",timeDiff, meanTime,vertex, Int_t(onlineMean), timeclock)); - } - //time in each channel as time[ipmt]-MeanTimeinThisChannel(with vertex=0) - for (Int_t ipmt=0; ipmt<24; ipmt++) { - if(time[ipmt]>1) { -// time[ipmt] = (time[ipmt] - fTime0vertex[ipmt])*channelWidth; - time[ipmt] = time[ipmt] * channelWidth; - frecpoints->SetTime(ipmt,time[ipmt]); - } - } + // } + clustersTree->Fill(); delete timeCFD; @@ -213,153 +199,157 @@ void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const { // T0 raw -> - // T0RecPoint writing - - //Q->T-> coefficients !!!! should be measured!!! - Int_t allData[110][50]; - - TArrayI * timeCFD = new TArrayI(24); - TArrayI * timeLED = new TArrayI(24); - TArrayI * chargeQT0 = new TArrayI(24); - TArrayI * chargeQT1 = new TArrayI(24); + // + // reference amplitude and time ref. point from reco param - TString option = GetOption(); - AliDebug(1,Form("Option: %s\n", option.Data())); + Int_t refAmp = GetRecoParam()->GetRefAmp(); + Int_t refPoint = GetRecoParam()->GetRefPoint(); + Int_t allData[110][5]; + + Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24]; + Int_t timeDiff=999999, meanTime=999999, timeclock=999999; + Float_t c = 0.0299792458; // cm/ps + Float_t vertex = 9999999; + Int_t onlineMean=0; for (Int_t i0=0; i0<105; i0++) { - for (Int_t j0=0; j0<50; j0++) allData[i0][j0]=0; - } - - - AliDebug(10," before read data "); - AliT0RawReader myrawreader(rawReader); - if (!myrawreader.Next()) - AliDebug(1,Form(" no raw data found!! %i", myrawreader.Next())); - - for (Int_t i=0; i<105; i++) { - for (Int_t iHit=0; iHit<50; iHit++) - { - allData[i][iHit] = myrawreader.GetData(i,iHit); - } + for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0; } - AliT0Calibrator *calib = new AliT0Calibrator(); - - // Int_t mV2Mip = param->GetmV2Mip(); - //mV2Mip = param->GetmV2Mip(); - Float_t channelWidth = fParam->GetChannelWidth() ; - - Int_t meanT0 = fParam->GetMeanT0(); - - for (Int_t in=0; in<24; in++) - { - timeLED->AddAt(allData[in+1][0],in); - timeCFD->AddAt(allData[in+25][0],in); - chargeQT1->AddAt(allData[in+57][0],in); - chargeQT0->AddAt(allData[in+80][0],in); - AliDebug(10, Form(" readed Raw %i %i %i %i %i", in, timeLED->At(in),timeCFD->At(in),chargeQT0->At(in),chargeQT1->At(in))); - } - + Float_t besttimeA=9999999; Float_t besttimeC=9999999; Int_t pmtBestA=99999; Int_t pmtBestC=99999; - Float_t timeDiff=9999999, meanTime=0; - - + Float_t meanVertex = fParam->GetMeanVertex(); + AliT0RecPoint* frecpoints= new AliT0RecPoint (); - - recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1); + recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1); + + AliDebug(10," before read data "); + AliT0RawReader myrawreader(rawReader); - Float_t time[24], adc[24]; - for (Int_t ipmt=0; ipmt<24; ipmt++) { - if(timeCFD->At(ipmt)>0 ){ - - if(option == "pdc"){ - Double_t qt0 = Double_t(chargeQT0->At(ipmt)); - Double_t qt1 = Double_t(chargeQT1->At(ipmt)); - if((qt1-qt0)>0) adc[ipmt] = TMath::Exp( Double_t (channelWidth*(qt1-qt0)/1000)); - // time[ipmt] = channelWidth * (calib-> WalkCorrection( ipmt,qt1 , timeCFD->At(ipmt) ) ) ; - time[ipmt] = calib-> WalkCorrection( ipmt,Int_t(qt1) , timeCFD->At(ipmt) ) ; - Double_t sl = (timeLED->At(ipmt) - time[ipmt])*channelWidth; - Double_t qt=((TGraph*)fAmpLEDrec.At(ipmt))->Eval(sl/1000.); - frecpoints->SetTime(ipmt,time[ipmt]); - frecpoints->SetAmp(ipmt,adc[ipmt]); - frecpoints->SetAmpLED(ipmt,qt); - AliDebug(1,Form(" QTC %f mv, QTC %f MIPS time in chann %f time %f ",adc[ipmt], adc[ipmt]/50.,time[ipmt], time[ipmt]*channelWidth)); - } - if(option == "cosmic") { - if(ipmt == 15) continue; //skip crashed PMT - Float_t qt0 = Float_t(chargeQT0->At(ipmt)); - Float_t qt1 = Float_t(chargeQT1->At(ipmt)); - if((qt0-qt1)>0) adc[ipmt] = qt0-qt1; - time[ipmt] = calib-> WalkCorrection( ipmt, Int_t(adc[ipmt]), timeCFD->At(ipmt) ) ; - Double_t sl = timeLED->At(ipmt) - time[ipmt]; - Double_t qt=((TGraph*)fAmpLEDrec.At(ipmt))->Eval(sl); - frecpoints->SetTime(ipmt,time[ipmt]); - frecpoints->SetAmp(ipmt,adc[ipmt]); - frecpoints->SetAmpLED(ipmt,qt); - AliDebug(10,Form(" QTC %i , time in chann %i led %i ", - Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( qt))); - } - - } - else { - time[ipmt] = 0; - adc[ipmt] = 0; - } - } - - for (Int_t ipmt=0; ipmt<12; ipmt++){ - if(time[ipmt] > 1 ) { - if(time[ipmt]GetType(); + + if (!myrawreader.Next()) + AliDebug(1,Form(" no raw data found!!")); + else + { + if(type == 7) { //only physics + for (Int_t i=0; i<105; i++) { + for (Int_t iHit=0; iHit<5; iHit++) + { + allData[i][iHit] = myrawreader.GetData(i,iHit); + } } - } - } - for ( Int_t ipmt=12; ipmt<24; ipmt++){ - if(time[ipmt] > 1) { - if(time[ipmt]SetTimeBestA(Int_t(besttimeA)); - if( besttimeC != 9999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC)); - AliDebug(1,Form(" besttimeA %f ps, besttimeC %f ps",besttimeA, besttimeC)); - Float_t c = 0.0299792; // cm/ps - Float_t vertex = 99999; - if(besttimeA <9999999 && besttimeC < 9999999 ){ - timeDiff =( besttimeC - besttimeA) *channelWidth; - if(option == "pdc") - meanTime = (meanT0 - (besttimeA + besttimeC)/2) * channelWidth; - if(option == "cosmic") meanTime = (besttimeA + besttimeC)/2; - vertex = c*(timeDiff)/2. + (fdZonA - fdZonC)/2; - AliDebug(1,Form(" timeDiff %f ps, meanTime %f ps, vertex %f cm",timeDiff, meanTime,vertex )); - frecpoints->SetVertex(vertex); - frecpoints->SetMeanTime(Int_t(meanTime)); - } - //time in each channel as time[ipmt]-MeanTimeinThisChannel(with vertex=0) - /* - for (Int_t ipmt=0; ipmt<24; ipmt++) { - if(time[ipmt]>1) { - time[ipmt] = (time[ipmt] - fTime0vertex[ipmt])*channelWidth; - frecpoints->SetTime(ipmt,time[ipmt]); - } - } - */ - recTree->Fill(); - - - delete timeCFD; - delete timeLED; - delete chargeQT0; - delete chargeQT1; - + Int_t ref = allData[refPoint][0]-5000.; + Float_t channelWidth = fParam->GetChannelWidth() ; + + // Int_t meanT0 = fParam->GetMeanT0(); + + + for (Int_t in=0; in<12; in++) + { + timeCFD[in] = allData[in+1][0] ; + timeCFD[in+12] = allData[in+56+1][0] ; + timeLED[in] = allData[in+12+1][0] ; + timeLED[in+12] = allData[in+68+1][0] ; + AliDebug(10, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ", + in, timeCFD[in],timeCFD[in+12],timeLED[in], + timeLED[in+12])); + } + + for (Int_t in=0; in<12; in++) + { + chargeQT0[in]=allData[2*in+25][0]; + chargeQT1[in]=allData[2*in+26][0]; + } + + for (Int_t in=12; in<24; in++) + { + chargeQT0[in]=allData[2*in+57][0]; + chargeQT1[in]=allData[2*in+58][0]; + } + + // } //cosmic with physics event + for (Int_t in=0; in<24; in++) + AliDebug(10, Form(" readed Raw %i %i %i %i %i", + in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in])); + onlineMean = allData[49][0]; + + Float_t time[24], adc[24]; + for (Int_t ipmt=0; ipmt<24; ipmt++) { + if(timeCFD[ipmt]>0 && timeLED[ipmt]>0){ + //for simulated data + //for physics data + if(( chargeQT1[ipmt] - chargeQT0[ipmt])>0) + adc[ipmt] = chargeQT1[ipmt] - chargeQT0[ipmt]; + else + adc[ipmt] = 0; + + + time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD[ipmt] ) ; + + Double_t sl = timeLED[ipmt] - timeCFD[ipmt]; + // time[ipmt] = fCalib-> WalkCorrection( ipmt, Int_t(sl), timeCFD[ipmt] ) ; + AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ", + ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl))); + Double_t ampMip =( (TGraph*)fAmpLED.At(ipmt))->Eval(sl); + Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]); + AliDebug(10,Form(" Amlitude in MIPS LED %f ; QTC %f; in channels %i\n ",ampMip,qtMip, adc[ipmt])); + + frecpoints->SetTime(ipmt, Float_t(time[ipmt]) ); + frecpoints->SetAmp(ipmt, Float_t( ampMip)); //for cosmic &pp beam + frecpoints->SetAmpLED(ipmt, Float_t(qtMip)); + + } + else { + time[ipmt] = 0; + adc[ipmt] = 0; + } + } + + for (Int_t ipmt=0; ipmt<12; ipmt++){ + if(time[ipmt] > 1 ) { + if(time[ipmt] 1) { + if(time[ipmt]SetTimeBestA(Int_t(besttimeA)); + if( besttimeC != 9999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC)); + AliDebug(5,Form(" pmtA %i besttimeA %f ps, pmtC %i besttimeC %f #channel", + pmtBestA,besttimeA, pmtBestC, besttimeC)); + if(besttimeA <9999999 && besttimeC < 9999999 ){ + timeDiff = ( besttimeA - besttimeC) *channelWidth; + meanTime = Float_t((besttimeA + besttimeC)/2.); + timeclock = Float_t(meanTime - ref); + vertex = meanVertex - c*(timeDiff)/2.; //+ (fdZonA - fdZonC)/2; + } + } //if phys event + frecpoints->SetT0clock(timeclock); + frecpoints->SetVertex(vertex); + frecpoints->SetMeanTime(Int_t(meanTime)); + frecpoints->SetOnlineMean(Int_t(onlineMean)); + AliDebug(5,Form(" timeDiff %f #channel, meanTime %f #channel, vertex %f cm meanVertex %f online mean %i ",timeDiff, meanTime,vertex,meanVertex, onlineMean)); + + + } // if (else )raw data + recTree->Fill(); + if(frecpoints) delete frecpoints; } + + //____________________________________________________________ void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const @@ -392,12 +382,14 @@ void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, Ali Float_t amp[24], time[24]; Float_t zPosition = frecpoints -> GetVertex(); Float_t timeStart = frecpoints -> GetMeanTime() ; - for ( Int_t i=0; i<24; i++) { + Float_t timeClock = frecpoints -> GetT0clock() ; + for ( Int_t i=0; i<24; i++) { time[i] = Float_t (frecpoints -> GetTime(i)); // ps to ns amp[i] = frecpoints -> GetAmp(i); } pESD->SetT0zVertex(zPosition); //vertex Z position pESD->SetT0(timeStart); // interaction time + pESD->SetT0clock(timeClock); // interaction time with ref.point(spectrum) pESD->SetT0time(time); // best TOF on each PMT pESD->SetT0amplitude(amp); // number of particles(MIPs) on each PMT