X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDCalibraFillHisto.cxx;h=4e0b3c68468ded28a3a6883c463c1cf7894b8b56;hb=132a16a18e5ed04406033e07e58309aa44ad5d39;hp=b79658b25cff16cc11c2a14840d0c048b1c3b16e;hpb=0c349049476aa6f7d12dc49f0502690879442e3e;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDCalibraFillHisto.cxx b/TRD/AliTRDCalibraFillHisto.cxx index b79658b25cf..4e0b3c68468 100644 --- a/TRD/AliTRDCalibraFillHisto.cxx +++ b/TRD/AliTRDCalibraFillHisto.cxx @@ -1,17 +1,17 @@ -/************************************************************************** - * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * * - * Author: The ALICE Off-line Project. * - * Contributors are mentioned in the code where appropriate. * - * * - * Permission to use, copy, modify and distribute this software and its * - * documentation strictly for non-commercial purposes is hereby granted * - * without fee, provided that the above copyright notice appears in all * - * copies and that both the copyright notice and this permission notice * - * appear in the supporting documentation. The authors make no claims * - * about the suitability of this software for any purpose. It is * - * provided "as is" without express or implied warranty. * - **************************************************************************/ +//************************************************************************** +// * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * +// * * +// * Author: The ALICE Off-line Project. * +// * Contributors are mentioned in the code where appropriate. * +// * * +// * Permission to use, copy, modify and distribute this software and its * +// * documentation strictly for non-commercial purposes is hereby granted * +// * without fee, provided that the above copyright notice appears in all * +// * copies and that both the copyright notice and this permission notice * +//* appear in the supporting documentation. The authors make no claims * +//* about the suitability of this software for any purpose. It is * +//* provided "as is" without express or implied warranty. * +//**************************************************************************/ /* $Id$ */ @@ -48,6 +48,7 @@ #include #include #include +#include #include "AliLog.h" @@ -57,11 +58,11 @@ #include "AliTRDCalibraVdriftLinearFit.h" #include "AliTRDcalibDB.h" #include "AliTRDCommonParam.h" -#include "AliTRDmcmTracklet.h" #include "AliTRDpadPlane.h" #include "AliTRDcluster.h" #include "AliTRDtrack.h" -#include "AliTRDRawStream.h" +#include "AliTRDtrackV1.h" +#include "AliTRDrawStreamBase.h" #include "AliRawReader.h" #include "AliRawReaderDate.h" #include "AliTRDgeometry.h" @@ -118,9 +119,7 @@ void AliTRDCalibraFillHisto::Terminate() AliTRDCalibraFillHisto::AliTRDCalibraFillHisto() :TObject() ,fGeo(0) - ,fMITracking(kFALSE) - ,fMcmTracking(kFALSE) - ,fMcmCorrectAngle(kFALSE) + ,fIsHLT(kFALSE) ,fCH2dOn(kFALSE) ,fPH2dOn(kFALSE) ,fPRF2dOn(kFALSE) @@ -130,31 +129,32 @@ AliTRDCalibraFillHisto::AliTRDCalibraFillHisto() ,fLinearFitterDebugOn(kFALSE) ,fRelativeScale(0) ,fThresholdClusterPRF2(15.0) + ,fLimitChargeIntegration(kFALSE) + ,fFillWithZero(kFALSE) + ,fNormalizeNbOfCluster(kFALSE) + ,fMaxCluster(0) + ,fNbMaxCluster(0) ,fCalibraMode(new AliTRDCalibraMode()) ,fDebugStreamer(0) ,fDebugLevel(0) - ,fDetectorAliTRDtrack(kFALSE) ,fDetectorPreviousTrack(-1) - ,fNumberClusters(18) + ,fMCMPrevious(-1) + ,fROBPrevious(-1) + ,fNumberClusters(1) + ,fNumberClustersf(30) ,fProcent(6.0) ,fDifference(17) ,fNumberTrack(0) ,fTimeMax(0) ,fSf(10.0) - ,fNumberBinCharge(100) - ,fNumberBinPRF(40) - ,fNgroupprf(0) - ,fListClusters(new TObjArray()) - ,fPar0(0x0) - ,fPar1(0x0) - ,fPar2(0x0) - ,fPar3(0x0) - ,fPar4(0x0) + ,fNumberBinCharge(50) + ,fNumberBinPRF(10) + ,fNgroupprf(3) ,fAmpTotal(0x0) ,fPHPlace(0x0) ,fPHValue(0x0) ,fGoodTracklet(kTRUE) - ,fGoodTrack(kTRUE) + ,fLinearFitterTracklet(0x0) ,fEntriesCH(0x0) ,fEntriesLinearFitter(0x0) ,fCalibraVector(0x0) @@ -165,8 +165,6 @@ AliTRDCalibraFillHisto::AliTRDCalibraFillHisto() ,fLinearVdriftFit(0x0) ,fCalDetGain(0x0) ,fCalROCGain(0x0) - ,fCalDetT0(0x0) - ,fCalROCT0(0x0) { // // Default constructor @@ -189,9 +187,7 @@ AliTRDCalibraFillHisto::AliTRDCalibraFillHisto() AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c) :TObject(c) ,fGeo(0) - ,fMITracking(c.fMITracking) - ,fMcmTracking(c.fMcmTracking) - ,fMcmCorrectAngle(c.fMcmCorrectAngle) + ,fIsHLT(c.fIsHLT) ,fCH2dOn(c.fCH2dOn) ,fPH2dOn(c.fPH2dOn) ,fPRF2dOn(c.fPRF2dOn) @@ -201,12 +197,19 @@ AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c) ,fLinearFitterDebugOn(c.fLinearFitterDebugOn) ,fRelativeScale(c.fRelativeScale) ,fThresholdClusterPRF2(c.fThresholdClusterPRF2) + ,fLimitChargeIntegration(c.fLimitChargeIntegration) + ,fFillWithZero(c.fFillWithZero) + ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster) + ,fMaxCluster(c.fMaxCluster) + ,fNbMaxCluster(c.fNbMaxCluster) ,fCalibraMode(0x0) ,fDebugStreamer(0) ,fDebugLevel(c.fDebugLevel) - ,fDetectorAliTRDtrack(c.fDetectorAliTRDtrack) ,fDetectorPreviousTrack(c.fDetectorPreviousTrack) + ,fMCMPrevious(c.fMCMPrevious) + ,fROBPrevious(c.fROBPrevious) ,fNumberClusters(c.fNumberClusters) + ,fNumberClustersf(c.fNumberClustersf) ,fProcent(c.fProcent) ,fDifference(c.fDifference) ,fNumberTrack(c.fNumberTrack) @@ -215,17 +218,11 @@ AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c) ,fNumberBinCharge(c.fNumberBinCharge) ,fNumberBinPRF(c.fNumberBinPRF) ,fNgroupprf(c.fNgroupprf) - ,fListClusters(new TObjArray()) - ,fPar0(0x0) - ,fPar1(0x0) - ,fPar2(0x0) - ,fPar3(0x0) - ,fPar4(0x0) ,fAmpTotal(0x0) ,fPHPlace(0x0) ,fPHValue(0x0) ,fGoodTracklet(c.fGoodTracklet) - ,fGoodTrack(c.fGoodTrack) + ,fLinearFitterTracklet(0x0) ,fEntriesCH(0x0) ,fEntriesLinearFitter(0x0) ,fCalibraVector(0x0) @@ -236,8 +233,6 @@ AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c) ,fLinearVdriftFit(0x0) ,fCalDetGain(0x0) ,fCalROCGain(0x0) - ,fCalDetT0(0x0) - ,fCalROCT0(0x0) { // // Copy constructor @@ -261,9 +256,7 @@ AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c) } if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain); - if(c.fCalDetT0) fCalDetT0 = new AliTRDCalDet(*c.fCalDetT0); if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain); - if(c.fCalROCT0) fCalROCT0 = new AliTRDCalROC(*c.fCalROCT0); if (fGeo) { delete fGeo; @@ -282,16 +275,25 @@ AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto() if ( fDebugStreamer ) delete fDebugStreamer; if ( fCalDetGain ) delete fCalDetGain; - if ( fCalDetT0 ) delete fCalDetT0; if ( fCalROCGain ) delete fCalROCGain; - if ( fCalROCT0 ) delete fCalROCT0; + if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; } + + delete [] fPHPlace; + delete [] fPHValue; + delete [] fEntriesCH; + delete [] fEntriesLinearFitter; + delete [] fAmpTotal; + + for(Int_t idet=0; idetGetNumberOfTimeBins(); - fSf = parCom->GetSamplingFrequency(); - fRelativeScale = 20; + fTimeMax = cal->GetNumberOfTimeBins(); + fSf = parCom->GetSamplingFrequency(); + if(!fNormalizeNbOfCluster) fRelativeScale = 20.0; + else fRelativeScale = 1.18; + fNumberClustersf = fTimeMax; + fNumberClusters = (Int_t)(0.5*fTimeMax); - //calib object from database used for reconstruction - if(fCalDetGain) delete fCalDetGain; - fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet())); - if(fCalDetT0) delete fCalDetT0; - fCalDetT0 = new AliTRDCalDet(*(cal->GetT0Det())); + // Init linear fitter + if(!fLinearFitterTracklet) { + fLinearFitterTracklet = new TLinearFitter(2,"pol1"); + fLinearFitterTracklet->StoreData(kTRUE); + } + //calib object from database used for reconstruction + if( fCalDetGain ){ + fCalDetGain->~AliTRDCalDet(); + new(fCalDetGain) AliTRDCalDet(*(cal->GetGainFactorDet())); + }else fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet())); + // Calcul Xbins Chambd0, Chamb2 Int_t ntotal0 = CalculateTotalNumberOfBins(0); Int_t ntotal1 = CalculateTotalNumberOfBins(1); @@ -411,23 +403,10 @@ Bool_t AliTRDCalibraFillHisto::Init2Dhistostrack() fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k)); fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k)); } - TString namech("Nz"); - namech += fCalibraMode->GetNz(0); - namech += "Nrphi"; - namech += fCalibraMode->GetNrphi(0); - fCalibraVector->SetNameCH((const char* ) namech); - TString nameph("Nz"); - nameph += fCalibraMode->GetNz(1); - nameph += "Nrphi"; - nameph += fCalibraMode->GetNrphi(1); - fCalibraVector->SetNamePH((const char* ) nameph); - TString nameprf("Nz"); - nameprf += fCalibraMode->GetNz(2); - nameprf += "Nrphi"; - nameprf += fCalibraMode->GetNrphi(2); - nameprf += "Ngp"; - nameprf += fNgroupprf; - fCalibraVector->SetNamePRF((const char* ) nameprf); + fCalibraVector->SetNzNrphi(0,fCalibraMode->GetNz(0),fCalibraMode->GetNrphi(0)); + fCalibraVector->SetNzNrphi(1,fCalibraMode->GetNz(1),fCalibraMode->GetNrphi(1)); + fCalibraVector->SetNzNrphi(2,fCalibraMode->GetNz(2),fCalibraMode->GetNrphi(2)); + fCalibraVector->SetNbGroupPRF(fNgroupprf); } // Create the 2D histos corresponding to the pad groupCalibration mode @@ -497,130 +476,47 @@ Bool_t AliTRDCalibraFillHisto::Init2Dhistostrack() return kTRUE; } - -//____________Functions for filling the histos in the code_____________________ - -//____________Offine tracking in the AliTRDtracker_____________________________ -Bool_t AliTRDCalibraFillHisto::ResetTrack() -{ - // - // For the offline tracking - // This function will be called in the function - // AliTRDtracker::FollowBackPropagation() at the beginning - // Reset the parameter to know we have a new TRD track - // - - fDetectorAliTRDtrack = kFALSE; - //fGoodTrack = kTRUE; - return kTRUE; - -} -//____________Offine tracking in the AliTRDtracker_____________________________ -void AliTRDCalibraFillHisto::ResetfVariables() -{ - // - // Reset values per tracklet - // - - // Reset the list of clusters - fListClusters->Clear(); - - //Reset the tracklet parameters - for(Int_t k = 0; k < fTimeMax; k++){ - fPar0[k] = 0.0; - fPar1[k] = 0.0; - fPar2[k] = 0.0; - fPar3[k] = 0.0; - fPar4[k] = 0.0; - } - - ResetfVariablestrack(); - -} -//____________Offine tracking in the AliTRDtracker_____________________________ -void AliTRDCalibraFillHisto::ResetfVariablestrack() -{ - // - // Reset values per tracklet - // - - //Reset good tracklet - fGoodTracklet = kTRUE; - - // Reset the fPHValue - if (fPH2dOn) { - //Reset the fPHValue and fPHPlace - for (Int_t k = 0; k < fTimeMax; k++) { - fPHValue[k] = 0.0; - fPHPlace[k] = -1; - } - } - - // Reset the fAmpTotal where we put value - if (fCH2dOn) { - for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) { - fAmpTotal[k] = 0.0; - } - } -} //____________Offline tracking in the AliTRDtracker____________________________ -Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t) +Bool_t AliTRDCalibraFillHisto::UpdateHistograms(const AliTRDtrack *t) { // - // For the offline tracking - // This function will be called in the function - // AliTRDtracker::FollowBackPropagation() in the loop over the clusters - // of TRD tracks - // Fill the 2D histos or the vectors with the info of the clusters at - // the end of a detectors if the track is "good" + // Use AliTRDtrack for the calibration // - - - AliTRDcluster *cl = 0x0; - Int_t index0 = 0; - Int_t index1 = 0; - // reset if good track - fGoodTrack = kTRUE; + AliTRDcluster *cl = 0x0; // pointeur to cluster + Int_t index0 = 0; // index of the first cluster in the current chamber + Int_t index1 = 0; // index of the current cluster in the current chamber - + ////////////////////////////// // loop over the clusters + /////////////////////////////// while((cl = t->GetCluster(index1))){ - // Localisation of the detector + ///////////////////////////////////////////////////////////////////////// + // Fill the infos for the previous clusters if not the same detector + //////////////////////////////////////////////////////////////////////// Int_t detector = cl->GetDetector(); - - - // Fill the infos for the previous clusters if not the same - // detector anymore but this time it should be the same track if ((detector != fDetectorPreviousTrack) && (index0 != index1)) { fNumberTrack++; - //printf("detector %d, fPreviousdetector %d, plane %d, planeprevious %d, index0 %d, index1 %d la\n",detector,fDetectorPreviousTrack,GetPlane(detector),GetPlane(fDetectorPreviousTrack),index0,index1); - //If the same track, then look if the previous detector is in //the same plane, if yes: not a good track - //FollowBack - //if (fDetectorAliTRDtrack && - // (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) { - //Follow - if ((GetPlane(detector) >= GetPlane(fDetectorPreviousTrack))) { - fGoodTrack = kFALSE; + if ((GetLayer(detector) == GetLayer(fDetectorPreviousTrack))) { + return kFALSE; } // Fill only if the track doesn't touch a masked pad or doesn't - // appear in the middle (fGoodTrack) - if (fGoodTrack && fGoodTracklet) { + if (fGoodTracklet) { // drift velocity unables to cut bad tracklets - Bool_t pass = FindP1TrackPHtrack(t,index0,index1); + Bool_t pass = FindP1TrackPHtracklet(t,index0,index1); // Gain calibration if (fCH2dOn) { - FillTheInfoOfTheTrackCH(); + FillTheInfoOfTheTrackCH(index1-index0); } // PH calibration @@ -628,18 +524,21 @@ Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t) FillTheInfoOfTheTrackPH(); } - if(pass && fPRF2dOn) HandlePRFtrack(t,index0,index1); + if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1); - } // if a good track + } // if a good tracklet // reset stuff - ResetfVariablestrack(); + ResetfVariablestracklet(); index0 = index1; } // Fill at the end the charge - - // Calcul the position of the detector and take the calib objects + + + ///////////////////////////////////////////////////////////// + // Localise and take the calib gain object for the detector + //////////////////////////////////////////////////////////// if (detector != fDetectorPreviousTrack) { //Localise the detector @@ -653,43 +552,48 @@ Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t) } // Get calib objects - if( fCalROCGain ) delete fCalROCGain; - fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector))); - if( fCalROCT0 ) delete fCalROCT0; - fCalROCT0 = new AliTRDCalROC(*(cal->GetT0ROC(detector))); + if( fCalROCGain ){ + if(!fIsHLT){ + fCalROCGain->~AliTRDCalROC(); + new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector))); + } + }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector))); } // Reset the detectbjobsor fDetectorPreviousTrack = detector; - // Store the info bis of the tracklet - Int_t *rowcol = CalculateRowCol(cl); - CheckGoodTracklet(detector,rowcol); + /////////////////////////////////////// + // Store the info of the cluster + /////////////////////////////////////// + Int_t row = cl->GetPadRow(); + Int_t col = cl->GetPadCol(); + CheckGoodTrackletV1(cl); Int_t group[2] = {0,0}; - if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,rowcol); - if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,rowcol); - StoreInfoCHPHtrack(cl,t,index1,group,rowcol); + if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col); + if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col); + StoreInfoCHPHtrack(cl,t->GetClusterdQdl(index1),group,row,col); index1++; } // while on clusters + /////////////////////////// // Fill the last plane + ////////////////////////// if( index0 != index1 ){ - - //printf("fPreviousdetector %d, planeprevious %d, index0 %d, index1 %d li\n",fDetectorPreviousTrack,GetPlane(fDetectorPreviousTrack),index0,index1); fNumberTrack++; - if (fGoodTrack && fGoodTracklet) { + if (fGoodTracklet) { // drift velocity unables to cut bad tracklets - Bool_t pass = FindP1TrackPHtrack(t,index0,index1); + Bool_t pass = FindP1TrackPHtracklet(t,index0,index1); // Gain calibration if (fCH2dOn) { - FillTheInfoOfTheTrackCH(); + FillTheInfoOfTheTrackCH(index1-index0); } // PH calibration @@ -697,1000 +601,1192 @@ Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t) FillTheInfoOfTheTrackPH(); } - if(pass && fPRF2dOn) HandlePRFtrack(t,index0,index1); + if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1); - } // if a good track + } // if a good tracklet } // reset stuff - ResetfVariablestrack(); + ResetfVariablestracklet(); return kTRUE; } //____________Offline tracking in the AliTRDtracker____________________________ -Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDcluster *cl, AliTRDtrack *t) +Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(const AliTRDtrackV1 *t) { // - // For the offline tracking - // This function will be called in the function - // AliTRDtracker::FollowBackPropagation() in the loop over the clusters - // of TRD tracks - // Fill the 2D histos or the vectors with the info of the clusters at - // the end of a detectors if the track is "good" + // Use AliTRDtrackV1 for the calibration // - // Localisation of the detector - Int_t detector = cl->GetDetector(); - - // Fill the infos for the previous clusters if not the same - // detector anymore or if not the same track - if (((detector != fDetectorPreviousTrack) || (!fDetectorAliTRDtrack)) && - (fDetectorPreviousTrack != -1)) { - - fNumberTrack++; + + const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane + AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet + Bool_t newtr = kTRUE; // new track + + // Get cal + AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); + if (!cal) { + AliInfo("Could not get calibDB"); + return kFALSE; + } + + /////////////////////////// + // loop over the tracklet + /////////////////////////// + for(Int_t itr = 0; itr < 6; itr++){ - // If the same track, then look if the previous detector is in - // the same plane, if yes: not a good track - //FollowBack - if (fDetectorAliTRDtrack && - (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) { - //Follow - //if (fDetectorAliTRDtrack && - // (GetPlane(detector) >= GetPlane(fDetectorPreviousTrack))) { - fGoodTrack = kFALSE; + if(!(tracklet = t->GetTracklet(itr))) continue; + if(!tracklet->IsOK()) continue; + fNumberTrack++; + ResetfVariablestracklet(); + + ////////////////////////////////////////// + // localisation of the tracklet and dqdl + ////////////////////////////////////////// + Int_t layer = tracklet->GetPlane(); + Int_t ic = 0; + while(!(cl = tracklet->GetClusters(ic++))) continue; + Int_t detector = cl->GetDetector(); + if (detector != fDetectorPreviousTrack) { + // if not a new track + if(!newtr){ + // don't use the rest of this track if in the same plane + if (layer == GetLayer(fDetectorPreviousTrack)) { + //printf("bad tracklet, same layer for detector %d\n",detector); + break; + } + } + //Localise the detector bin + LocalisationDetectorXbins(detector); + // Get calib objects + if( fCalROCGain ){ + if(!fIsHLT){ + fCalROCGain->~AliTRDCalROC(); + new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector))); + } + }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector))); + + // reset + fDetectorPreviousTrack = detector; } - - // Fill only if the track doesn't touch a masked pad or doesn't - // appear in the middle (fGoodTrack) - if (fGoodTrack && fGoodTracklet) { + newtr = kFALSE; + + //////////////////////////// + // loop over the clusters + //////////////////////////// + Int_t nbclusters = 0; + for(int jc=0; jcGetClusters(jc))) continue; + nbclusters++; + + // Store the info bis of the tracklet + Int_t row = cl->GetPadRow(); + Int_t col = cl->GetPadCol(); + CheckGoodTrackletV1(cl); + Int_t group[2] = {0,0}; + if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col); + if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col); + StoreInfoCHPHtrack(cl, tracklet->GetdQdl(jc),group,row,col); + } + + //////////////////////////////////////// + // Fill the stuffs if a good tracklet + //////////////////////////////////////// + if (fGoodTracklet) { // drift velocity unables to cut bad tracklets - Bool_t pass = FindP1TrackPH(); - + Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters); + + //printf("pass %d and nbclusters %d\n",pass,nbclusters); + // Gain calibration if (fCH2dOn) { - FillTheInfoOfTheTrackCH(); + FillTheInfoOfTheTrackCH(nbclusters); } - + // PH calibration if (fPH2dOn) { FillTheInfoOfTheTrackPH(); } - - if(pass && fPRF2dOn) HandlePRF(); - - - } // if a good track - - ResetfVariables(); - if(!fDetectorAliTRDtrack) fGoodTrack = kTRUE; - - } // Fill at the end the charge - - // Calcul the position of the detector and take the calib objects - if (detector != fDetectorPreviousTrack) { - //Localise the detector - LocalisationDetectorXbins(detector); - - // Get cal - AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); - if (!cal) { - AliInfo("Could not get calibDB"); - return kFALSE; - } - - // Get calib objects - if( fCalROCGain ) delete fCalROCGain; - fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector))); - if( fCalROCT0 ) delete fCalROCT0; - fCalROCT0 = new AliTRDCalROC(*(cal->GetT0ROC(detector))); + + if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters); + + } // if a good tracklet } - - // Reset the detector - fDetectorPreviousTrack = detector; - fDetectorAliTRDtrack = kTRUE; - - // Store the infos of the tracklets - AliTRDcluster *kcl = new AliTRDcluster(*cl); - fListClusters->Add((TObject *)kcl); - Int_t time = cl->GetLocalTimeBin(); - fPar0[time] = t->GetY(); - fPar1[time] = t->GetZ(); - fPar2[time] = t->GetSnp(); - fPar3[time] = t->GetTgl(); - fPar4[time] = t->GetSigned1Pt(); - - // Store the info bis of the tracklet - Int_t *rowcol = CalculateRowCol(cl); - CheckGoodTracklet(detector,rowcol); - Int_t group[2] = {0,0}; - if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,rowcol); - if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,rowcol); - StoreInfoCHPH(cl,t,group,rowcol); - + return kTRUE; } -//____________Online trackling in AliTRDtrigger________________________________ -Bool_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk) +/////////////////////////////////////////////////////////////////////////////////// +// Routine inside the update with AliTRDtrack +/////////////////////////////////////////////////////////////////////////////////// +//____________Offine tracking in the AliTRDtracker_____________________________ +Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtracklet(const AliTRDtrack *t, Int_t index0, Int_t index1) { // - // For the tracking - // This function will be called in the function AliTRDtrigger::TestTracklet - // before applying the pt cut on the tracklets - // Fill the infos for the tracklets fTrkTest if the tracklets is "good" + // Drift velocity calibration: + // Fit the clusters with a straight line + // From the slope find the drift velocity // + + //Number of points: if less than 3 return kFALSE + Int_t npoints = index1-index0; + if(npoints <= 2) return kFALSE; + + ///////////////// + //Variables + //////////////// + Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector + // parameters of the track + Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in the plan yx track + Double_t tgl = t->GetTglPlane(GetLayer(detector)); // dz/dl and not dz/dx! + Double_t tnp = 0.0; // tan angle in the plan xy track + if( TMath::Abs(snp) < 1.){ + tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp)); + } + Float_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx now from dz/dl + // tilting pad and cross row + Int_t crossrow = 0; // if it crosses a pad row + Int_t rowp = -1; // if it crosses a pad row + AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(detector),GetStack(detector)); + Double_t tiltingangle = padplane->GetTiltingAngle(); // tiltingangle of the pad + Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle + // linear fit + fLinearFitterTracklet->ClearPoints(); + Double_t dydt = 0.0; // dydt tracklet after straight line fit + Double_t errorpar = 0.0; // error after straight line fit on dy/dt + Double_t pointError = 0.0; // error after straight line fit + Int_t nbli = 0; // number linear fitter points - // Localisation of the Xbins involved - Int_t idect = trk->GetDetector(); - fDetectorPreviousTrack = idect; - LocalisationDetectorXbins(idect); + ////////////////////////////// + // loop over clusters + //////////////////////////// + for(Int_t k = 0; k < npoints; k++){ + + AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0); + if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue; + Double_t ycluster = cl->GetY(); + Int_t time = cl->GetPadTime(); + Double_t timeis = time/fSf; + //Double_t q = cl->GetQ(); + //See if cross two pad rows + Int_t row = cl->GetPadRow(); + if(k==0) rowp = row; + if(row != rowp) crossrow = 1; + + fLinearFitterTracklet->AddPoint(&timeis,ycluster,1); + nbli++; - // Get the parameter object - AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); - if (!cal) { - AliInfo("Could not get calibDB"); - return kFALSE; } - - // Reset - ResetfVariables(); - - // Get calib objects - if( fCalROCGain ) delete fCalROCGain; - fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(idect))); - if( fCalROCT0 ) delete fCalROCT0; - fCalROCT0 = new AliTRDCalROC(*(cal->GetT0ROC(idect))); - - // Row of the tracklet and position in the pad groups - Int_t *rowcol = new Int_t[2]; - rowcol[0] = trk->GetRow(); - Int_t group[3] = {-1,-1,-1}; - - // Eventuelle correction due to track angle in z direction - Float_t correction = 1.0; - if (fMcmCorrectAngle) { - Float_t z = trk->GetRowz(); - Float_t r = trk->GetTime0(); - correction = r / TMath::Sqrt((r*r+z*z)); + + ////////////////////////////// + // linear fit + ////////////////////////////// + if(nbli <= 2){ + fLinearFitterTracklet->ClearPoints(); + return kFALSE; } + TVectorD pars; + fLinearFitterTracklet->Eval(); + fLinearFitterTracklet->GetParameters(pars); + pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2)); + errorpar = fLinearFitterTracklet->GetParError(1)*pointError; + dydt = pars[1]; + fLinearFitterTracklet->ClearPoints(); + + ///////////////////////////// + // debug + //////////////////////////// + if(fDebugLevel > 0){ + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + + (* fDebugStreamer) << "FindP1TrackPHtracklet0"<< + //"snpright="<GetNclusters() >= fNumberClusters) { - for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) { + Int_t nbclusters = index1-index0; + Int_t layer = GetLayer(fDetectorPreviousTrack); - Float_t amp[3] = { 0.0, 0.0, 0.0 }; - Int_t time = trk->GetClusterTime(icl); - rowcol[1] = trk->GetClusterCol(icl); - - amp[0] = trk->GetClusterADC(icl)[0] * correction; - amp[1] = trk->GetClusterADC(icl)[1] * correction; - amp[2] = trk->GetClusterADC(icl)[2] * correction; + (* fDebugStreamer) << "FindP1TrackPHtracklet1"<< + //"snpright="< fNumberClustersf) return kFALSE; + if(pointError >= 0.1) return kFALSE; + if(crossrow == 1) return kFALSE; - // Col of cluster and position in the pad groups - if(fCH2dOn) { - group[0] = CalculateCalibrationGroup(0,rowcol); - fAmpTotal[(Int_t) group[0]] += (Float_t) (amp[0]+amp[1]+amp[2]); - } - if(fPH2dOn) { - group[1] = CalculateCalibrationGroup(1,rowcol); - fPHPlace[time] = group[1]; - fPHValue[time] = (Float_t) (amp[0]+amp[1]+amp[2]); + //////////////////////////// + // fill + //////////////////////////// + if(fLinearFitterOn){ + //Add to the linear fitter of the detector + if( TMath::Abs(snp) < 1.){ + Double_t x = tnp-dzdx*tnt; + (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt); + if(fLinearFitterDebugOn) { + fLinearVdriftFit->Update(detector,x,pars[1]); } - if(fPRF2dOn) group[2] = CalculateCalibrationGroup(2,rowcol); - - // See if we are not near a masked pad fGoodTracklet - CheckGoodTracklet(idect,rowcol); - - // Fill PRF direct without tnp bins...only for monitoring... - if (fPRF2dOn && fGoodTracklet) { - - if ((amp[0] > fThresholdClusterPRF2) && - (amp[1] > fThresholdClusterPRF2) && - (amp[2] > fThresholdClusterPRF2) && - ((amp[0]*amp[2]/(amp[1]*amp[1])) < 0.06)) { - - // Security of the denomiateur is 0 - if ((((Float_t) (((Float_t) amp[1]) * ((Float_t) amp[1]))) - / ((Float_t) (((Float_t) amp[0]) * ((Float_t) amp[2])))) != 1.0) { - Float_t xcenter = 0.5 * (TMath::Log(amp[2] / amp[0])) - / (TMath::Log((amp[1]*amp[1]) / (amp[0]*amp[2]))); - Float_t ycenter = amp[1] / (amp[0] + amp[1] + amp[2]); - - if (TMath::Abs(xcenter) < 0.5) { - Float_t yminus = amp[0] / (amp[0]+amp[1]+amp[2]); - Float_t ymax = amp[2] / (amp[0]+amp[1]+amp[2]); - // Fill only if it is in the drift region! - if (((Float_t) time / fSf) > 0.3) { - if (fHisto2d) { - fPRF2d->Fill(xcenter,(fCalibraMode->GetXbins(2)+group[2]+0.5),ycenter); - fPRF2d->Fill(-(xcenter+1.0),(fCalibraMode->GetXbins(2)+group[2]+0.5),yminus); - fPRF2d->Fill((1.0-xcenter),(fCalibraMode->GetXbins(2)+group[2]+0.5),ymax); - } - if (fVector2d) { - fCalibraVector->UpdateVectorPRF(idect,group[2],xcenter,ycenter); - fCalibraVector->UpdateVectorPRF(idect,group[2],-(xcenter+1.0),yminus); - fCalibraVector->UpdateVectorPRF(idect,group[2],(1.0-xcenter),ymax); - } - }//in the drift region - }//in the middle - }//denominateur security - }//cluster shape and thresholds - }//good and PRF On - - } // Boucle clusters - - // Fill the charge - if(fGoodTracklet){ - if (fCH2dOn) FillTheInfoOfTheTrackCH(); - if (fPH2dOn) FillTheInfoOfTheTrackPH(); + fEntriesLinearFitter[detector]++; } - - fNumberTrack++; - - } // Condition on number of clusters - - return kTRUE; + } + return kTRUE; } -//_____________________________________________________________________________ -Int_t *AliTRDCalibraFillHisto::CalculateRowCol(AliTRDcluster *cl) const +//____________Offine tracking in the AliTRDtracker_____________________________ +Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters) { // - // Calculate the row and col number of the cluster + // Drift velocity calibration: + // Fit the clusters with a straight line + // From the slope find the drift velocity // + //////////////////////////////////////////////// + //Number of points: if less than 3 return kFALSE + ///////////////////////////////////////////////// + if(nbclusters <= 2) return kFALSE; - Int_t *rowcol = new Int_t[2]; - rowcol[0] = 0; - rowcol[1] = 0; - - // Localisation of the detector - Int_t detector = cl->GetDetector(); - Int_t chamber = GetChamber(detector); - Int_t plane = GetPlane(detector); - - // Localisation of the cluster - Double_t pos[3] = { 0.0, 0.0, 0.0 }; - pos[0] = ((AliCluster *)cl)->GetX(); - pos[1] = cl->GetY(); - pos[2] = cl->GetZ(); + //////////// + //Variables + //////////// + // results of the linear fit + Double_t dydt = 0.0; // dydt tracklet after straight line fit + Double_t errorpar = 0.0; // error after straight line fit on dy/dt + Double_t pointError = 0.0; // error after straight line fit + // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant + Int_t crossrow = 0; // if it crosses a pad row + Int_t rowp = -1; // if it crosses a pad row + Float_t tnt = tracklet->GetTilt(); // tan tiltingangle + fLinearFitterTracklet->ClearPoints(); + + + /////////////////////////////////////////// + // Take the parameters of the track + ////////////////////////////////////////// + // take now the snp, tnp and tgl from the track + Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber + Double_t tnp = 0.0; // dy/dx at the end of the chamber + if( TMath::Abs(snp) < 1.){ + tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp)); + } + Double_t tgl = tracklet->GetTgl(); // dz/dl + Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl + // at the entrance + //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber + //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber + //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl + // at the end with correction due to linear fit + //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction + //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction + + + //////////////////////////// + // loop over the clusters + //////////////////////////// + Int_t nbli = 0; + AliTRDcluster *cl = 0x0; + for(int ic=0; icGetClusters(ic))) continue; + if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue; + + Double_t ycluster = cl->GetY(); + Int_t time = cl->GetPadTime(); + Double_t timeis = time/fSf; + //See if cross two pad rows + Int_t row = cl->GetPadRow(); + if(rowp==-1) rowp = row; + if(row != rowp) crossrow = 1; - // Position of the cluster - AliTRDpadPlane *padplane = fGeo->GetPadPlane(plane,chamber); - Int_t row = padplane->GetPadRowNumber(pos[2]); - //Do not take from here because it was corrected from ExB already.... - //Double_t offsetz = padplane->GetPadRowOffset(row,pos[2]); - //Double_t offsettilt = padplane->GetTiltOffset(offsetz); - //Int_t col = padplane->GetPadColNumber(pos[1] + offsettilt,offsetz); - //Int_t col = padplane->GetPadColNumber(pos[1]+offsettilt); - Int_t col = cl->GetPadCol(); + fLinearFitterTracklet->AddPoint(&timeis,ycluster,1); + nbli++; - //return - rowcol[0] = row; - rowcol[1] = col; - return rowcol; + ////////////////////////////// + // Check no shared clusters + ////////////////////////////// + for(int icc=AliTRDseedV1::kNtb; iccGetClusters(icc))) crossrow = 1; + } + } -} -//_____________________________________________________________________________ -void AliTRDCalibraFillHisto::CheckGoodTracklet(Int_t detector, Int_t *rowcol) -{ - // - // See if we are not near a masked pad - // + //////////////////////////////////// + // Do the straight line fit now + /////////////////////////////////// + if(nbli <= 2){ + fLinearFitterTracklet->ClearPoints(); + return kFALSE; + } + TVectorD pars; + fLinearFitterTracklet->Eval(); + fLinearFitterTracklet->GetParameters(pars); + pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2)); + errorpar = fLinearFitterTracklet->GetParError(1)*pointError; + dydt = pars[1]; + //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar); + fLinearFitterTracklet->ClearPoints(); + + //////////////////////////////// + // Debug stuff + /////////////////////////////// - Int_t row = rowcol[0]; - Int_t col = rowcol[1]; - if (!IsPadOn(detector, col, row)) { - fGoodTracklet = kFALSE; - } - - if (col > 0) { - if (!IsPadOn(detector, col-1, row)) { - fGoodTracklet = kFALSE; - } - } + if(fDebugLevel > 0){ + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + - if (col < 143) { - if (!IsPadOn(detector, col+1, row)) { - fGoodTracklet = kFALSE; - } - } - -} -//_____________________________________________________________________________ -Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t col, Int_t row) const -{ - // - // Look in the choosen database if the pad is On. - // If no the track will be "not good" - // + Int_t layer = GetLayer(fDetectorPreviousTrack); + + (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<< + //"snpright="<IsChamberInstalled(detector) || - cal->IsChamberMasked(detector) || - cal->IsPadMasked(detector,col,row)) { - return kFALSE; - } - else { - return kTRUE; - } - -} -//_____________________________________________________________________________ -Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t *rowcol) const -{ - // - // Calculate the calibration group number for i - // - - // Row of the cluster and position in the pad groups - Int_t posr = 0; - if (fCalibraMode->GetNnZ(i) != 0) { - posr = (Int_t) rowcol[0] / fCalibraMode->GetNnZ(i); - } - - - // Col of the cluster and position in the pad groups - Int_t posc = 0; - if (fCalibraMode->GetNnRphi(i) != 0) { - posc = (Int_t) rowcol[1] / fCalibraMode->GetNnRphi(i); } - return posc*fCalibraMode->GetNfragZ(i)+posr; - -} -//____________________________________________________________________________________ -Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i) -{ - // - // Calculate the total number of calibration groups - // - - Int_t ntotal = 0; - fCalibraMode->ModePadCalibration(2,i); - fCalibraMode->ModePadFragmentation(0,2,0,i); - fCalibraMode->SetDetChamb2(i); - ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i); - fCalibraMode->ModePadCalibration(0,i); - fCalibraMode->ModePadFragmentation(0,0,0,i); - fCalibraMode->SetDetChamb0(i); - ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i); - AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i)); - return ntotal; + ///////////////////////// + // Cuts quality + //////////////////////// -} -//____________Set the pad calibration variables for the detector_______________ -Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector) -{ - // - // For the detector calcul the first Xbins and set the number of row - // and col pads per calibration groups, the number of calibration - // groups in the detector. - // + if(nbclusters < fNumberClusters) return kFALSE; + if(nbclusters > fNumberClustersf) return kFALSE; + if(pointError >= 0.3) return kFALSE; + if(crossrow == 1) return kFALSE; - // first Xbins of the detector - if (fCH2dOn) { - fCalibraMode->CalculXBins(detector,0); - } - if (fPH2dOn) { - fCalibraMode->CalculXBins(detector,1); - } - if (fPRF2dOn) { - fCalibraMode->CalculXBins(detector,2); - } + /////////////////////// + // Fill + ////////////////////// - // fragmentation of idect - for (Int_t i = 0; i < 3; i++) { - fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i); - fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector) - , (Int_t) GetChamber(detector) - , (Int_t) GetSector(detector),i); + if(fLinearFitterOn){ + //Add to the linear fitter of the detector + if( TMath::Abs(snp) < 1.){ + Double_t x = tnp-dzdx*tnt; + (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt); + if(fLinearFitterDebugOn) { + fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]); + } + fEntriesLinearFitter[fDetectorPreviousTrack]++; + } } return kTRUE; - } -//_____________________________________________________________________________ -void AliTRDCalibraFillHisto::StoreInfoCHPH(AliTRDcluster *cl, AliTRDtrack *t, Int_t *group, Int_t *rowcol) +//____________Offine tracking in the AliTRDtracker_____________________________ +Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(const AliTRDtrack *t, Int_t index0, Int_t index1) { // - // Store the infos in fAmpTotal, fPHPlace and fPHValue - // - - // Charge in the cluster - Float_t q = TMath::Abs(cl->GetQ()); - Int_t time = cl->GetLocalTimeBin(); - - //Correct for the gain coefficient used in the database for reconstruction - Float_t correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(rowcol[1],rowcol[0]); - Float_t correcttheT0 = fCalDetT0->GetValue(fDetectorPreviousTrack)+fCalROCT0->GetValue(rowcol[1],rowcol[0]); - - // we substract correcttheT0 in AliTRDclusterizerV1::MakeClusters (line 458) - Int_t timec = Arrondi((Double_t)(time+correcttheT0)); - if((correcttheT0+0.5)==(int(correcttheT0+0.5))) { - timec++; - } - if( timec < 0 ) return; - - - // Correction due to the track angle - Float_t correction = 1.0; - Float_t normalisation = 6.67; - // we divide with gain in AliTRDclusterizerV1::Transform... - if( correctthegain > 0 ) normalisation /= correctthegain; - if ((q >0) && (t->GetNdedx() > 0)) { - correction = t->GetClusterdQdl((t->GetNdedx() - 1)) / (normalisation); - } - - // Fill the fAmpTotal with the charge - if (fCH2dOn) { - fAmpTotal[(Int_t) group[0]] += correction; - } - - // Fill the fPHPlace and value - if (fPH2dOn) { - fPHPlace[timec] = group[1]; - fPHValue[timec] = correction; - } - -} -//_____________________________________________________________________________ -void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, AliTRDtrack *t, Int_t index, Int_t *group, Int_t *rowcol) -{ + // PRF width calibration + // Assume a Gaussian shape: determinate the position of the three pad clusters + // Fit with a straight line + // Take the fitted values for all the clusters (3 or 2 pad clusters) + // Fill the PRF as function of angle of the track // - // Store the infos in fAmpTotal, fPHPlace and fPHValue // - - // Charge in the cluster - Float_t q = TMath::Abs(cl->GetQ()); - Int_t time = cl->GetLocalTimeBin(); - - //Correct for the gain coefficient used in the database for reconstruction - Float_t correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(rowcol[1],rowcol[0]); - Float_t correcttheT0 = fCalDetT0->GetValue(fDetectorPreviousTrack)+fCalROCT0->GetValue(rowcol[1],rowcol[0]); - - // we substract correcttheT0 in AliTRDclusterizerV1::MakeClusters (line 458) - Int_t timec = Arrondi((Double_t)(time+correcttheT0)); - if((correcttheT0+0.5)==(int(correcttheT0+0.5))) { - timec++; - } - if( timec < 0 ) return; - // Correction due to the track angle - Float_t correction = 1.0; - Float_t normalisation = 6.67; - // we divide with gain in AliTRDclusterizerV1::Transform... - if( correctthegain > 0 ) normalisation /= correctthegain; - if (q >0) { - correction = t->GetClusterdQdl(index) / (normalisation); - } - // Fill the fAmpTotal with the charge - if (fCH2dOn) { - fAmpTotal[(Int_t) group[0]] += correction; + ////////////////////////// + // variables + ///////////////////////// + Int_t npoints = index1-index0; // number of total points + Int_t nb3pc = 0; // number of three pads clusters used for fit + Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector + // To see the difference due to the fit + Double_t *padPositions; + padPositions = new Double_t[npoints]; + for(Int_t k = 0; k < npoints; k++){ + padPositions[k] = 0.0; + } + // Take the tgl and snp with the track t now + Double_t tgl = t->GetTglPlane(GetLayer(detector)); //dz/dl and not dz/dx + Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in xy plan + Float_t dzdx = 0.0; // dzdx + Float_t tnp = 0.0; + if(TMath::Abs(snp) < 1.0){ + tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp)); + dzdx = tgl*TMath::Sqrt(1+tnp*tnp); } - - // Fill the fPHPlace and value - if (fPH2dOn) { - fPHPlace[timec] = group[1]; - fPHValue[timec] = correction; + // linear fitter + fLinearFitterTracklet->ClearPoints(); + + /////////////////////////// + // calcul the tnp group + /////////////////////////// + Bool_t echec = kFALSE; + Double_t shift = 0.0; + //Calculate the shift in x coresponding to this tnp + if(fNgroupprf != 0.0){ + shift = -3.0*(fNgroupprf-1)-1.5; + Double_t limithigh = -0.2*(fNgroupprf-1); + if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE; + else{ + while(tnp > limithigh){ + limithigh += 0.2; + shift += 3.0; + } + } } - -} -//_____________________________________________________________________ -Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDRawStream *rawStream, Bool_t nocheck) -{ - // - // Event Processing loop - AliTRDRawStream - // 0 timebin problem - // 1 no input - // 2 input - // - - Int_t withInput = 1; - - Int_t phvalue[36]; - //Int_t row[36]; - //Int_t col[36]; - for(Int_t k = 0; k < 36; k++){ - phvalue[k] = 10; - //row[k] = -1; - //col[36] = -1; + if(echec) { + delete [] padPositions; + return kFALSE; } - fDetectorPreviousTrack = -1; - Int_t nbtimebin = 0; - if(!nocheck){ - - fTimeMax = 0; - - while (rawStream->Next()) { - - Int_t idetector = rawStream->GetDet(); // current detector - if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){ - if(TMath::Mean(fTimeMax,phvalue)>20.0){ - withInput = 2; - for(Int_t k = 0; k < fTimeMax; k++){ - UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],fTimeMax); - phvalue[k] = 10; - //row[k] = -1; - //col[k] = -1; - } - } - } - fDetectorPreviousTrack = idetector; - nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data - if(nbtimebin == 0) return 0; - if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0; - fTimeMax = nbtimebin; - Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin - //row[iTimeBin] = rawStream->GetRow(); // current row - //col[iTimeBin] = rawStream->GetCol(); // current col - Int_t *signal = rawStream->GetSignals(); // current ADC signal - - Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3)); - Int_t n = 0; - for(Int_t itime = iTimeBin; itime < fin; itime++){ - // should extract baseline here! - if(signal[n]>13) phvalue[itime] = signal[n]; - n++; + ////////////////////// + // loop clusters + ///////////////////// + for(Int_t k = 0; k < npoints; k++){ + //Take the cluster + AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0); + if(!cl) continue; + Short_t *signals = cl->GetSignals(); + Double_t time = cl->GetPadTime(); + //Calculate x if possible + Float_t xcenter = 0.0; + Bool_t echec1 = kTRUE; + if((time<=7) || (time>=21)) continue; + // Center 3 balanced: position with the center of the pad + if ((((Float_t) signals[3]) > 0.0) && + (((Float_t) signals[2]) > 0.0) && + (((Float_t) signals[4]) > 0.0)) { + echec1 = kFALSE; + // Security if the denomiateur is 0 + if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) / + ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) { + xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2])))) + / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) + / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4]))))); } - } - - // fill the last one - if(fDetectorPreviousTrack != -1){ - if(TMath::Mean(fTimeMax,phvalue)>20.0){ - withInput = 2; - for(Int_t k = 0; k < fTimeMax; k++){ - UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],fTimeMax); - phvalue[k] = 10; - //row[k] = -1; - //col[k] = -1; - } + else { + echec1 = kTRUE; } } - - } - else{ + if(TMath::Abs(xcenter) > 0.5) echec = kTRUE; + if(echec) continue; + //if no echec: calculate with the position of the pad + // Position of the cluster + Double_t padPosition = xcenter + cl->GetPadCol(); + padPositions[k] = padPosition; + nb3pc++; + fLinearFitterTracklet->AddPoint(&time, padPosition,1); + }//clusters loop - while (rawStream->Next()) { - Int_t idetector = rawStream->GetDet(); // current detector - if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){ - if(TMath::Mean(nbtimebin,phvalue)>20.0){ - withInput = 2; - for(Int_t k = 0; k < nbtimebin; k++){ - UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin); - phvalue[k] = 10; - //row[k] = -1; - //col[k] = -1; - } - } - } - fDetectorPreviousTrack = idetector; - nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data - Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin - //row[iTimeBin] = rawStream->GetRow(); // current row - //col[iTimeBin] = rawStream->GetCol(); // current col - Int_t *signal = rawStream->GetSignals(); // current ADC signal + ///////////////////////////// + // fit + //////////////////////////// + if(nb3pc < 3) { + delete [] padPositions; + fLinearFitterTracklet->ClearPoints(); + return kFALSE; + } + fLinearFitterTracklet->Eval(); + TVectorD line(2); + fLinearFitterTracklet->GetParameters(line); + Float_t pointError = -1.0; + if( fLinearFitterTracklet->GetChisquare()>=0.0) { + pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2)); + } + fLinearFitterTracklet->ClearPoints(); + + ///////////////////////////////////////////////////// + // Now fill the PRF: second loop over clusters + ///////////////////////////////////////////////////// + for(Int_t k = 0; k < npoints; k++){ + //Take the cluster + AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0); + Short_t *signals = cl->GetSignals(); // signal + Double_t time = cl->GetPadTime(); // time bin + Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit + Float_t padPos = cl->GetPadCol(); // middle pad + Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit + Float_t ycenter = 0.0; // relative center charge + Float_t ymin = 0.0; // relative left charge + Float_t ymax = 0.0; // relative right charge + + //Requiere simply two pads clusters at least + if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) || + ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){ + Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]); + if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum; + if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum; + if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum; + } + + //calibration group + Int_t rowcl = cl->GetPadRow(); // row of cluster + Int_t colcl = cl->GetPadCol(); // col of cluster + Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group + Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group + Float_t xcl = cl->GetY(); // y cluster + Float_t qcl = cl->GetQ(); // charge cluster + Int_t layer = GetLayer(detector); // layer + Int_t stack = GetStack(detector); // stack + Double_t xdiff = dpad; // reconstructed position constant + Double_t x = dpad; // reconstructed position moved + Float_t ep = pointError; // error of fit + Float_t signal1 = (Float_t)signals[1]; // signal at the border + Float_t signal3 = (Float_t)signals[3]; // signal + Float_t signal2 = (Float_t)signals[2]; // signal + Float_t signal4 = (Float_t)signals[4]; // signal + Float_t signal5 = (Float_t)signals[5]; // signal at the border + + ////////////////////////////// + // debug + ///////////////////////////// + if(fDebugLevel > 0){ + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + + x = xdiff; + Int_t type=0; + Float_t y = ycenter; + (* fDebugStreamer) << "HandlePRFtracklet"<< + "caligroup="<13) phvalue[itime] = signal[n]; - n++; + } + + //////////////////////////// + // quality cuts + /////////////////////////// + if(npoints < fNumberClusters) continue; + if(npoints > fNumberClustersf) continue; + if(nb3pc <= 5) continue; + if((time >= 21) || (time < 7)) continue; + if(TMath::Abs(snp) >= 1.0) continue; + if(TMath::Abs(qcl) < 80) continue; + + //////////////////////////// + // Fill + /////////////////////////// + if (fHisto2d) { + if(TMath::Abs(dpad) < 1.5) { + fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter); + fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter); + } + if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) { + fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin); + fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin); + } + if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) { + fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax); + fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax); } } - - // fill the last one - if(fDetectorPreviousTrack != -1){ - if(TMath::Mean(nbtimebin,phvalue)>20.0){ - withInput = 2; - for(Int_t k = 0; k < nbtimebin; k++){ - UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin); - phvalue[k] = 10; - //row[k] = -1; - //col[k] = -1; - } + if (fVector2d) { + if(TMath::Abs(dpad) < 1.5) { + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter); + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter); + } + if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) { + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin); + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin); + } + if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) { + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax); + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax); } } } + delete [] padPositions; + return kTRUE; - return withInput; - } -//_____________________________________________________________________ -Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck) +//____________Offine tracking in the AliTRDtracker_____________________________ +Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters) { // - // Event processing loop - AliRawReader + // PRF width calibration + // Assume a Gaussian shape: determinate the position of the three pad clusters + // Fit with a straight line + // Take the fitted values for all the clusters (3 or 2 pad clusters) + // Fill the PRF as function of angle of the track + // // + //printf("begin\n"); + /////////////////////////////////////////// + // Take the parameters of the track + ////////////////////////////////////////// + // take now the snp, tnp and tgl from the track + Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber + Double_t tnp = 0.0; // dy/dx at the end of the chamber + if( TMath::Abs(snp) < 1.){ + tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp)); + } + Double_t tgl = tracklet->GetTgl(); // dz/dl + Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl + // at the entrance + //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber + //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber + //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl + // at the end with correction due to linear fit + //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction + //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction + + /////////////////////////////// + // Calculate tnp group shift + /////////////////////////////// + Bool_t echec = kFALSE; + Double_t shift = 0.0; + //Calculate the shift in x coresponding to this tnp + if(fNgroupprf != 0.0){ + shift = -3.0*(fNgroupprf-1)-1.5; + Double_t limithigh = -0.2*(fNgroupprf-1); + if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE; + else{ + while(tnp > limithigh){ + limithigh += 0.2; + shift += 3.0; + } + } + } + // do nothing if out of tnp range + //printf("echec %d\n",(Int_t)echec); + if(echec) return kFALSE; - AliTRDRawStream rawStream(rawReader); + /////////////////////// + // Variables + ////////////////////// - rawReader->Select("TRD"); + Int_t nb3pc = 0; // number of three pads clusters used for fit + // to see the difference between the fit and the 3 pad clusters position + Double_t padPositions[AliTRDseedV1::kNtb]; + memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t)); + fLinearFitterTracklet->ClearPoints(); + + //printf("loop clusters \n"); + //////////////////////////// + // loop over the clusters + //////////////////////////// + AliTRDcluster *cl = 0x0; + for(int ic=0; icGetClusters(ic+AliTRDseedV1::kNtb))) continue; + } + if(!(cl = tracklet->GetClusters(ic))) continue; + Double_t time = cl->GetPadTime(); + if((time<=7) || (time>=21)) continue; + Short_t *signals = cl->GetSignals(); + Float_t xcenter = 0.0; + Bool_t echec1 = kTRUE; - return ProcessEventDAQ(&rawStream, nocheck); -} -//_________________________________________________________________________ -Int_t AliTRDCalibraFillHisto::ProcessEventDAQ( -#ifdef ALI_DATE - eventHeaderStruct *event, - Bool_t nocheck -#else - eventHeaderStruct* /*event*/, - Bool_t /*nocheck*/ - -#endif - ) -{ - // - // process date event - // -#ifdef ALI_DATE - AliRawReader *rawReader = new AliRawReaderDate((void*)event); - Int_t result=ProcessEventDAQ(rawReader, nocheck); - delete rawReader; - return result; -#else - Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE"); - return 0; -#endif + ///////////////////////////////////////////////////////////// + // Center 3 balanced: position with the center of the pad + ///////////////////////////////////////////////////////////// + if ((((Float_t) signals[3]) > 0.0) && + (((Float_t) signals[2]) > 0.0) && + (((Float_t) signals[4]) > 0.0)) { + echec1 = kFALSE; + // Security if the denomiateur is 0 + if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) / + ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) { + xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2])))) + / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) + / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4]))))); + } + else { + echec1 = kTRUE; + } + } + if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE; + if(echec1) continue; -} -//____________Online trackling in AliTRDtrigger________________________________ -Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Int_t signal, Int_t nbtimebins) -{ - // - // For the DAQ - // Fill a simple average pulse height - // - - // Localisation of the Xbins involved - //LocalisationDetectorXbins(det); + //////////////////////////////////////////////////////// + //if no echec1: calculate with the position of the pad + // Position of the cluster + // fill the linear fitter + /////////////////////////////////////////////////////// + Double_t padPosition = xcenter + cl->GetPadCol(); + padPositions[ic] = padPosition; + nb3pc++; + fLinearFitterTracklet->AddPoint(&time, padPosition,1); - // Row and position in the pad groups - //Int_t posr = 0; - //if (fCalibraMode->GetNnZ(1) != 0) { - // posr = (Int_t) row / fCalibraMode->GetNnZ(1); - //} + + }//clusters loop + + //printf("Fin loop clusters \n"); + ////////////////////////////// + // fit with a straight line + ///////////////////////////// + if(nb3pc < 3){ + fLinearFitterTracklet->ClearPoints(); + return kFALSE; + } + fLinearFitterTracklet->Eval(); + TVectorD line(2); + fLinearFitterTracklet->GetParameters(line); + Float_t pointError = -1.0; + if( fLinearFitterTracklet->GetChisquare()>=0.0) { + pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2)); + } + fLinearFitterTracklet->ClearPoints(); - // Col of cluster and position in the pad groups - //Int_t posc = 0; - //if (fCalibraMode->GetNnRphi(1) != 0) { - // posc = (Int_t) col / fCalibraMode->GetNnRphi(1); - //} - - //fPH2d->Fill((Float_t) timebin/fSf,(fCalibraMode->GetXbins(1)+posc*fCalibraMode->GetNfragZ(1)+posr)+0.5,(Float_t) signal); - - ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal); - - return kTRUE; + //printf("PRF second loop \n"); + //////////////////////////////////////////////// + // Fill the PRF: Second loop over clusters + ////////////////////////////////////////////// + for(int ic=0; icGetClusters(ic+AliTRDseedV1::kNtb))) continue; + // + if(!(cl = tracklet->GetClusters(ic))) continue; + + Short_t *signals = cl->GetSignals(); // signal + Double_t time = cl->GetPadTime(); // time bin + Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit + Float_t padPos = cl->GetPadCol(); // middle pad + Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit + Float_t ycenter = 0.0; // relative center charge + Float_t ymin = 0.0; // relative left charge + Float_t ymax = 0.0; // relative right charge -} -//____________Write_____________________________________________________ -//_____________________________________________________________________ -void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append) + //////////////////////////////////////////////////////////////// + // Calculate ycenter, ymin and ymax even for two pad clusters + //////////////////////////////////////////////////////////////// + if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) || + ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){ + Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]); + if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum; + if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum; + if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum; + } + + ///////////////////////// + // Calibration group + //////////////////////// + Int_t rowcl = cl->GetPadRow(); // row of cluster + Int_t colcl = cl->GetPadCol(); // col of cluster + Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group + Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group + Float_t xcl = cl->GetY(); // y cluster + Float_t qcl = cl->GetQ(); // charge cluster + Int_t layer = GetLayer(fDetectorPreviousTrack); // layer + Int_t stack = GetStack(fDetectorPreviousTrack); // stack + Double_t xdiff = dpad; // reconstructed position constant + Double_t x = dpad; // reconstructed position moved + Float_t ep = pointError; // error of fit + Float_t signal1 = (Float_t)signals[1]; // signal at the border + Float_t signal3 = (Float_t)signals[3]; // signal + Float_t signal2 = (Float_t)signals[2]; // signal + Float_t signal4 = (Float_t)signals[4]; // signal + Float_t signal5 = (Float_t)signals[5]; // signal at the border + + + + ///////////////////// + // Debug stuff + //////////////////// + + if(fDebugLevel > 0){ + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + x = xdiff; + Int_t type=0; + Float_t y = ycenter; + (* fDebugStreamer) << "HandlePRFtrackletV1"<< + "caligroup="< fNumberClustersf) continue; + if(nb3pc <= 5) continue; + if((time >= 21) || (time < 7)) continue; + if(TMath::Abs(qcl) < 80) continue; + if( TMath::Abs(snp) > 1.) continue; + + + //////////////////////// + // Fill the histos + /////////////////////// + if (fHisto2d) { + if(TMath::Abs(dpad) < 1.5) { + fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter); + fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter); + //printf("place %f, ycenter %f\n",(shift+dpad),ycenter); + } + if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) { + fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin); + fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin); + } + if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) { + fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax); + fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax); + } + } + // vector method + if (fVector2d) { + if(TMath::Abs(dpad) < 1.5) { + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter); + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter); + } + if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) { + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin); + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin); + } + if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) { + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax); + fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax); + } + } + } // second loop over clusters + + + return kTRUE; +} +/////////////////////////////////////////////////////////////////////////////////////// +// Pad row col stuff: see if masked or not +/////////////////////////////////////////////////////////////////////////////////////// +//_____________________________________________________________________________ +void AliTRDCalibraFillHisto::CheckGoodTrackletV1(const AliTRDcluster *cl) { // - // Write infos to file + // See if we are not near a masked pad // - - //For debugging - if ( fDebugStreamer ) { - delete fDebugStreamer; - fDebugStreamer = 0x0; - } - AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d" - ,fNumberTrack - ,fNumberUsedCh[0] - ,fNumberUsedCh[1] - ,fNumberUsedPh[0] - ,fNumberUsedPh[1])); - - TDirectory *backup = gDirectory; - TString option; - - if ( append ) - option = "update"; - else - option = "recreate"; - - TFile f(filename,option.Data()); + if(cl->IsMasked()) fGoodTracklet = kFALSE; + - TStopwatch stopwatch; - stopwatch.Start(); - if(fVector2d) { - f.WriteTObject(fCalibraVector); - } +} +//_____________________________________________________________________________ +void AliTRDCalibraFillHisto::CheckGoodTrackletV0(const Int_t detector,const Int_t row,const Int_t col) +{ + // + // See if we are not near a masked pad + // - if (fCH2dOn ) { - if (fHisto2d) { - f.WriteTObject(fCH2d); - } + if (!IsPadOn(detector, col, row)) { + fGoodTracklet = kFALSE; } - if (fPH2dOn ) { - if (fHisto2d) { - f.WriteTObject(fPH2d); + + if (col > 0) { + if (!IsPadOn(detector, col-1, row)) { + fGoodTracklet = kFALSE; } } - if (fPRF2dOn) { - if (fHisto2d) { - f.WriteTObject(fPRF2d); + + if (col < 143) { + if (!IsPadOn(detector, col+1, row)) { + fGoodTracklet = kFALSE; } } - if(fLinearFitterOn){ - AnalyseLinearFitter(); - f.WriteTObject(fLinearVdriftFit); + +} +//_____________________________________________________________________________ +Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const +{ + // + // Look in the choosen database if the pad is On. + // If no the track will be "not good" + // + + // Get the parameter object + AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); + if (!cal) { + AliInfo("Could not get calibDB"); + return kFALSE; } - - f.Close(); - if ( backup ) backup->cd(); + if (!cal->IsChamberInstalled(detector) || + cal->IsChamberMasked(detector) || + cal->IsPadMasked(detector,col,row)) { + return kFALSE; + } + else { + return kTRUE; + } - AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs" - ,stopwatch.RealTime(),stopwatch.CpuTime())); } -//___________________________________________probe the histos__________________________________________________ -Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i) +/////////////////////////////////////////////////////////////////////////////////////// +// Calibration groups: calculate the number of groups, localise... +//////////////////////////////////////////////////////////////////////////////////////// +//_____________________________________________________________________________ +Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const { // - // Check the number of stats in h, 0 is TH2I 1 is TProfile2D - // debug mode with 2 for TH2I and 3 for TProfile2D - // It gives a pointer to a Double_t[7] with the info following... - // [0] : number of calibration groups with entries - // [1] : minimal number of entries found - // [2] : calibration group number of the min - // [3] : maximal number of entries found - // [4] : calibration group number of the max - // [5] : mean number of entries found - // [6] : mean relative error + // Calculate the calibration group number for i // - - Double_t *info = new Double_t[7]; - - // Number of Xbins (detectors or groups of pads) - Int_t nbins = h->GetNbinsY(); //number of calibration groups - Int_t nxbins = h->GetNbinsX(); //number of bins per histo - - // Initialise - Double_t nbwe = 0; //number of calibration groups with entries - Double_t minentries = 0; //minimal number of entries found - Double_t maxentries = 0; //maximal number of entries found - Double_t placemin = 0; //calibration group number of the min - Double_t placemax = -1; //calibration group number of the max - Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry - Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D - - Double_t counter = 0; - - //Debug - TH1F *nbEntries = 0x0;//distribution of the number of entries - TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group - TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule - - // Beginning of the loop over the calibration groups - for (Int_t idect = 0; idect < nbins; idect++) { - - TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e"); - projch->SetDirectory(0); - - // Number of entries for this calibration group - Double_t nentries = 0.0; - if((i%2) == 0){ - for (Int_t k = 0; k < nxbins; k++) { - nentries += h->GetBinContent(h->GetBin(k+1,idect+1)); - } - } - else{ - for (Int_t k = 0; k < nxbins; k++) { - nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1)); - if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) { - meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1))))); - counter++; - } - } - } - - //Debug - if(i > 1){ - if((!((Bool_t)nbEntries)) && (nentries > 0)){ - nbEntries = new TH1F("Number of entries","Number of entries" - ,100,(Int_t)nentries/2,nentries*2); - nbEntries->SetDirectory(0); - nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group" - ,nbins,0,nbins); - nbEntriesPerGroup->SetDirectory(0); - nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule" - ,(Int_t)(nbins/18),0,(Int_t)(nbins/18)); - nbEntriesPerSp->SetDirectory(0); - } - if(nbEntries){ - if(nentries > 0) nbEntries->Fill(nentries); - nbEntriesPerGroup->Fill(idect+0.5,nentries); - nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries); - } - } - - //min amd max - if(nentries > maxentries){ - maxentries = nentries; - placemax = idect; - } - if(idect == 0) { - minentries = nentries; - } - if(nentries < minentries){ - minentries = nentries; - placemin = idect; - } - //nbwe - if(nentries > 0) { - nbwe++; - meanstats += nentries; - } - }//calibration groups loop - - if(nbwe > 0) meanstats /= nbwe; - if(counter > 0) meanrelativerror /= counter; - - AliInfo(Form("There are %f calibration groups with entries",nbwe)); - AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin)); - AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax)); - AliInfo(Form("The mean number of entries is %f",meanstats)); - if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror)); - - info[0] = nbwe; - info[1] = minentries; - info[2] = placemin; - info[3] = maxentries; - info[4] = placemax; - info[5] = meanstats; - info[6] = meanrelativerror; - - if(i > 1){ - gStyle->SetPalette(1); - gStyle->SetOptStat(1111); - gStyle->SetPadBorderMode(0); - gStyle->SetCanvasColor(10); - gStyle->SetPadLeftMargin(0.13); - gStyle->SetPadRightMargin(0.01); - TCanvas *stat = new TCanvas("stat","",50,50,600,800); - stat->Divide(2,1); - stat->cd(1); - nbEntries->Draw(""); - stat->cd(2); - nbEntriesPerSp->SetStats(0); - nbEntriesPerSp->Draw(""); - TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800); - stat1->cd(); - nbEntriesPerGroup->SetStats(0); - nbEntriesPerGroup->Draw(""); + + // Row of the cluster and position in the pad groups + Int_t posr = 0; + if (fCalibraMode->GetNnZ(i) != 0) { + posr = (Int_t) row / fCalibraMode->GetNnZ(i); } - - return info; - -} -//____________________________________________________________________________ -Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH() -{ - // - // Return a Int_t[4] with: - // 0 Mean number of entries - // 1 median of number of entries - // 2 rms of number of entries - // 3 number of group with entries - // - - Double_t *stat = new Double_t[4]; - stat[3] = 0.0; - - Int_t nbofgroups = CalculateTotalNumberOfBins(0); - Double_t *weight = new Double_t[nbofgroups]; - Int_t *nonul = new Int_t[nbofgroups]; - for(Int_t k = 0; k < nbofgroups; k++){ - if(fEntriesCH[k] > 0) { - weight[k] = 1.0; - nonul[(Int_t)stat[3]] = fEntriesCH[k]; - stat[3]++; - } - else weight[k] = 0.0; + + // Col of the cluster and position in the pad groups + Int_t posc = 0; + if (fCalibraMode->GetNnRphi(i) != 0) { + posc = (Int_t) col / fCalibraMode->GetNnRphi(i); } - stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight); - stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight); - stat[2] = TMath::RMS((Int_t)stat[3],nonul); - - return stat; - + + return posc*fCalibraMode->GetNfragZ(i)+posr; + } -//____________________________________________________________________________ -Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const +//____________________________________________________________________________________ +Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i) { // - // Return a Int_t[4] with: - // 0 Mean number of entries - // 1 median of number of entries - // 2 rms of number of entries - // 3 number of group with entries + // Calculate the total number of calibration groups // + + Int_t ntotal = 0; - Double_t *stat = new Double_t[4]; - stat[3] = 0.0; - - Int_t nbofgroups = 540; - Double_t *weight = new Double_t[nbofgroups]; - Int_t *nonul = new Int_t[nbofgroups]; - - for(Int_t k = 0; k < nbofgroups; k++){ - if(fEntriesLinearFitter[k] > 0) { - weight[k] = 1.0; - nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k]; - stat[3]++; - } - else weight[k] = 0.0; - } - stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight); - stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight); - stat[2] = TMath::RMS((Int_t)stat[3],nonul); - - return stat; - -} -//_____________________________________________________________________________ -void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF) -{ - // - // Should be between 0 and 6 - // - - if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) { - AliInfo("The number of groups must be between 0 and 6!"); - } - else { - fNgroupprf = numberGroupsPRF; + // All together + if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){ + ntotal = 1; + AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i)); + return ntotal; } -} -//_____________________________________________________________________________ -void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale) -{ - // - // Set the factor that will divide the deposited charge - // to fit in the histo range [0,300] - // - - if (RelativeScale > 0.0) { - fRelativeScale = RelativeScale; - } - else { - AliInfo("RelativeScale must be strict positif!"); + // Per Supermodule + if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){ + ntotal = 18; + AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i)); + return ntotal; } -} + // More + fCalibraMode->ModePadCalibration(2,i); + fCalibraMode->ModePadFragmentation(0,2,0,i); + fCalibraMode->SetDetChamb2(i); + ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i); + fCalibraMode->ModePadCalibration(0,i); + fCalibraMode->ModePadFragmentation(0,0,0,i); + fCalibraMode->SetDetChamb0(i); + ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i); + AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i)); + return ntotal; +} //_____________________________________________________________________________ void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz) { @@ -1707,7 +1803,6 @@ void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz) } } - //_____________________________________________________________________________ void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi) { @@ -1724,89 +1819,155 @@ void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi) } } -//____________Protected Functions______________________________________________ -//____________Create the 2D histo to be filled online__________________________ -// + //_____________________________________________________________________________ -void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn) +void AliTRDCalibraFillHisto::SetAllTogether(Int_t i) { // - // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each - // If fNgroupprf is zero then no binning in tnp + // Set the mode of calibration group all together // - - TString name("Nz"); - name += fCalibraMode->GetNz(2); - name += "Nrphi"; - name += fCalibraMode->GetNrphi(2); - name += "Ngp"; - name += fNgroupprf; - - if(fNgroupprf != 0){ - - fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name - ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn); - fPRF2d->SetYTitle("Det/pad groups"); - fPRF2d->SetXTitle("Position x/W [pad width units]"); - fPRF2d->SetZTitle("Q_{i}/Q_{total}"); - fPRF2d->SetStats(0); + if(fVector2d == kTRUE) { + AliInfo("Can not work with the vector method"); + return; } - else{ - fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name - ,fNumberBinPRF,-1.5,1.5,nn,0,nn); - fPRF2d->SetYTitle("Det/pad groups"); - fPRF2d->SetXTitle("Position x/W [pad width units]"); - fPRF2d->SetZTitle("Q_{i}/Q_{total}"); - fPRF2d->SetStats(0); - } - + fCalibraMode->SetAllTogether(i); + } //_____________________________________________________________________________ -void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn) +void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i) { // - // Create the 2D histos + // Set the mode of calibration group per supermodule // + if(fVector2d == kTRUE) { + AliInfo("Can not work with the vector method"); + return; + } + fCalibraMode->SetPerSuperModule(i); + +} - TString name("Nz"); - name += fCalibraMode->GetNz(1); - name += "Nrphi"; - name += fCalibraMode->GetNrphi(1); +//____________Set the pad calibration variables for the detector_______________ +Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector) +{ + // + // For the detector calcul the first Xbins and set the number of row + // and col pads per calibration groups, the number of calibration + // groups in the detector. + // - fPH2d = new TProfile2D("PH2d",(const Char_t *) name - ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf - ,nn,0,nn); - fPH2d->SetYTitle("Det/pad groups"); - fPH2d->SetXTitle("time [#mus]"); - fPH2d->SetZTitle(" [a.u.]"); - fPH2d->SetStats(0); + // first Xbins of the detector + if (fCH2dOn) { + fCalibraMode->CalculXBins(detector,0); + } + if (fPH2dOn) { + fCalibraMode->CalculXBins(detector,1); + } + if (fPRF2dOn) { + fCalibraMode->CalculXBins(detector,2); + } + + // fragmentation of idect + for (Int_t i = 0; i < 3; i++) { + fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i); + fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector) + , (Int_t) GetStack(detector) + , (Int_t) GetSector(detector),i); + } + + return kTRUE; } //_____________________________________________________________________________ -void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn) +void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF) { // - // Create the 2D histos + // Should be between 0 and 6 // + + if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) { + AliInfo("The number of groups must be between 0 and 6!"); + } + else { + fNgroupprf = numberGroupsPRF; + } - TString name("Nz"); - name += fCalibraMode->GetNz(0); - name += "Nrphi"; - name += fCalibraMode->GetNrphi(0); +} +/////////////////////////////////////////////////////////////////////////////////////////// +// Per tracklet: store or reset the info, fill the histos with the info +////////////////////////////////////////////////////////////////////////////////////////// +//_____________________________________________________________________________ +void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(const AliTRDcluster *cl,const Double_t dqdl,const Int_t *group,const Int_t row,const Int_t col) +{ + // + // Store the infos in fAmpTotal, fPHPlace and fPHValue + // Correct from the gain correction before + // + + // time bin of the cluster not corrected + Int_t time = cl->GetPadTime(); + Float_t charge = TMath::Abs(cl->GetQ()); - fCH2d = new TH2I("CH2d",(const Char_t *) name - ,fNumberBinCharge,0,300,nn,0,nn); - fCH2d->SetYTitle("Det/pad groups"); - fCH2d->SetXTitle("charge deposit [a.u]"); - fCH2d->SetZTitle("counts"); - fCH2d->SetStats(0); - fCH2d->Sumw2(); + //printf("Store::time %d, amplitude %f\n",time,dqdl); + + //Correct for the gain coefficient used in the database for reconstruction + Float_t correctthegain = 1.0; + if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack); + else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row); + Float_t correction = 1.0; + Float_t normalisation = 6.67; + // we divide with gain in AliTRDclusterizer::Transform... + if( correctthegain > 0 ) normalisation /= correctthegain; + + + // take dd/dl corrected from the angle of the track + correction = dqdl / (normalisation); + + // Fill the fAmpTotal with the charge + if (fCH2dOn) { + if((!fLimitChargeIntegration) || (cl->IsInChamber())) { + //printf("Store::group %d, amplitude %f\n",group[0],correction); + fAmpTotal[(Int_t) group[0]] += correction; + } + } + + // Fill the fPHPlace and value + if (fPH2dOn) { + fPHPlace[time] = group[1]; + fPHValue[time] = charge; + } + } +//____________Offine tracking in the AliTRDtracker_____________________________ +void AliTRDCalibraFillHisto::ResetfVariablestracklet() +{ + // + // Reset values per tracklet + // + + //Reset good tracklet + fGoodTracklet = kTRUE; + + // Reset the fPHValue + if (fPH2dOn) { + //Reset the fPHValue and fPHPlace + for (Int_t k = 0; k < fTimeMax; k++) { + fPHValue[k] = 0.0; + fPHPlace[k] = -1; + } + } + // Reset the fAmpTotal where we put value + if (fCH2dOn) { + for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) { + fAmpTotal[k] = 0.0; + } + } +} //____________Offine tracking in the AliTRDtracker_____________________________ -void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH() +void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters) { // // For the offline tracking or mcm tracklets @@ -1818,11 +1979,21 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH() Int_t fd = -1; // Premiere zone non nulle Float_t totalcharge = 0.0; // Total charge for the supermodule histo - - + //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf); + + if(nbclusters < fNumberClusters) return; + if(nbclusters > fNumberClustersf) return; + + + // Normalize with the number of clusters + Double_t normalizeCst = fRelativeScale; + if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters; + + //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)); // See if the track goes through different zones for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) { + //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k); if (fAmpTotal[k] > 0.0) { totalcharge += fAmpTotal[k]; nb++; @@ -1832,18 +2003,19 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH() } } - + //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack); + switch (nb) { case 1: fNumberUsedCh[0]++; fEntriesCH[fCalibraMode->GetXbins(0)+fd]++; if (fHisto2d) { - FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale); - //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5); + FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst); + //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5); } if (fVector2d) { - fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale); + fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst); } break; case 2: @@ -1852,22 +2024,22 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH() // One of the two very big if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) { if (fHisto2d) { - FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale); - //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5); + FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst); + //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5); } if (fVector2d) { - fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale); + fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst); } fNumberUsedCh[1]++; fEntriesCH[fCalibraMode->GetXbins(0)+fd]++; } if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) { if (fHisto2d) { - FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale); - //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5); + FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst); + //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5); } if (fVector2d) { - fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/fRelativeScale); + fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst); } fNumberUsedCh[1]++; fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++; @@ -1880,24 +2052,24 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH() // One of the two very big if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) { if (fHisto2d) { - FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale); - //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5); + FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst); + //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5); } if (fVector2d) { - fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale); + fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst); } fNumberUsedCh[1]++; fEntriesCH[fCalibraMode->GetXbins(0)+fd]++; } if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) { if (fHisto2d) { - FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale); - //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5); + FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst); + //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5); } fNumberUsedCh[1]++; fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++; if (fVector2d) { - fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale); + fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst); } } } @@ -1922,10 +2094,14 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH() Int_t fd2 = -1; // Deuxieme zone non nulle Int_t k1 = -1; // Debut de la premiere zone Int_t k2 = -1; // Debut de la seconde zone + Int_t nbclusters = 0; // number of clusters + + // See if the track goes through different zones for (Int_t k = 0; k < fTimeMax; k++) { if (fPHValue[k] > 0.0) { + nbclusters++; if (fd1 == -1) { fd1 = fPHPlace[k]; k1 = k; @@ -1943,6 +2119,20 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH() } } + // See if noise before and after + if(fMaxCluster > 0) { + if(fPHValue[0] > fMaxCluster) return; + if(fTimeMax > fNbMaxCluster) { + for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){ + if(fPHValue[k] > fMaxCluster) return; + } + } + } + + //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf); + + if(nbclusters < fNumberClusters) return; + if(nbclusters > fNumberClustersf) return; switch(nb) { @@ -1950,10 +2140,17 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH() fNumberUsedPh[0]++; for (Int_t i = 0; i < fTimeMax; i++) { if (fHisto2d) { - fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + } + //printf("Fill the time bin %d with %f\n",i,fPHValue[i]); } if (fVector2d) { - fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + } } } break; @@ -1966,10 +2163,16 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH() fNumberUsedPh[1]++; for (Int_t i = 0; i < fTimeMax; i++) { if (fHisto2d) { - fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + } } if (fVector2d) { - fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + } } } } @@ -1978,10 +2181,16 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH() fNumberUsedPh[1]++; for (Int_t i = 0; i < fTimeMax; i++) { if (fHisto2d) { - fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]); + if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]); + } } if (fVector2d) { - fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]); + if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]); + } } } } @@ -1997,10 +2206,16 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH() fNumberUsedPh[1]++; for (Int_t i = 0; i < fTimeMax; i++) { if (fHisto2d) { - fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + } } if (fVector2d) { - fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + } } } } @@ -2009,10 +2224,16 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH() fNumberUsedPh[1]++; for (Int_t i = 0; i < fTimeMax; i++) { if (fHisto2d) { - fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]); + if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]); + } } if (fVector2d) { - fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]); + if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]); + } } } } @@ -2028,10 +2249,16 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH() fNumberUsedPh[1]++; for (Int_t i = 0; i < fTimeMax; i++) { if (fHisto2d) { - fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]); + } } if (fVector2d) { - fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]); + } } } } @@ -2040,10 +2267,16 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH() fNumberUsedPh[1]++; for (Int_t i = 0; i < fTimeMax; i++) { if (fHisto2d) { - fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]); + if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]); + } } if (fVector2d) { - fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]); + if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]); + else { + if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]); + } } } } @@ -2054,904 +2287,792 @@ void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH() default: break; } } -//____________Offine tracking in the AliTRDtracker_____________________________ -Bool_t AliTRDCalibraFillHisto::FindP1TrackPH() +////////////////////////////////////////////////////////////////////////////////////////// +// DAQ process functions +///////////////////////////////////////////////////////////////////////////////////////// +//_____________________________________________________________________ +Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamBase *rawStream, Bool_t nocheck) { // - // For the offline tracking - // This function will be called in the functions UpdateHistogram... - // to fill the find the parameter P1 of a track for the drift velocity calibration + // Event Processing loop - AliTRDrawStreamBase + // TestBeam 2007 version + // 0 timebin problem + // 1 no input + // 2 input + // Same algorithm as TestBeam but different reader // - - - //Number of points: if less than 3 return kFALSE - Int_t npoints = fListClusters->GetEntriesFast(); - if(npoints <= 2) return kFALSE; - - //Variables - TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet - Double_t snp = 0.0; // sin angle in the plan yx track - Double_t y = 0.0; // y clusters in the middle of the chamber - Double_t z = 0.0; // z cluster in the middle of the chamber - Double_t dydt = 0.0; // dydt tracklet after straight line fit - Double_t tnp = 0.0; // tan angle in the plan xy track - Double_t tgl = 0.0; // dz/dl and not dz/dx! - Double_t errorpar = 0.0; // error after straight line fit on dy/dt - Double_t pointError = 0.0; // error after straight line fit - Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); //detector - Int_t snpright = 1; // if we took in the middle snp - Int_t crossrow = 0; // if it crosses a pad row - Double_t tiltingangle = 0; // tiltingangle of the pad - Float_t dzdx = 0; // dz/dx now from dz/dl - Int_t nbli = 0; // number linear fitter points - AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector)); - - linearFitterTracklet.StoreData(kFALSE); - linearFitterTracklet.ClearPoints(); - //if more than one row - Int_t rowp = -1; // if it crosses a pad row - - //tiltingangle - tiltingangle = padplane->GetTiltingAngle(); - Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle - - //Fill with points - for(Int_t k = 0; k < npoints; k++){ - - AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k); - Double_t ycluster = cl->GetY(); - Int_t time = cl->GetLocalTimeBin(); - Double_t timeis = time/fSf; - //See if cross two pad rows - Int_t row = padplane->GetPadRowNumber(cl->GetZ()); - if(k==0) rowp = row; - if(row != rowp) crossrow = 1; - //Take in the middle of the chamber - //FollowBack - if(time > (Int_t) 10) { - //Follow - //if(time < (Int_t) 11) { - z = cl->GetZ(); - y = cl->GetY(); - snp = fPar2[time]; - tgl = fPar3[time]; + Int_t withInput = 1; + + Double_t phvalue[16][144][36]; + for(Int_t k = 0; k < 36; k++){ + for(Int_t j = 0; j < 16; j++){ + for(Int_t c = 0; c < 144; c++){ + phvalue[j][c][k] = 0.0; + } } - linearFitterTracklet.AddPoint(&timeis,ycluster,1); - nbli++; } - //FollowBack - if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() < 10) snpright = 0; - //Follow - //if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() >= 11) snpright = 0; - if(nbli <= 2) return kFALSE; - // Do the straight line fit now - TVectorD pars; - linearFitterTracklet.Eval(); - linearFitterTracklet.GetParameters(pars); - pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli); - errorpar = linearFitterTracklet.GetParError(1)*pointError; - dydt = pars[1]; - - if( TMath::Abs(snp) < 1.){ - tnp = snp / (TMath::Sqrt(1-(snp*snp))); - } - dzdx = tgl*TMath::Sqrt(1+tnp*tnp); + fDetectorPreviousTrack = -1; + fMCMPrevious = -1; + fROBPrevious = -1; + Int_t nbtimebin = 0; + Int_t baseline = 10; - if(fDebugLevel > 0){ - if ( !fDebugStreamer ) { - //debug stream - TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root"); - if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer - } - - (* fDebugStreamer) << "VDRIFT0"<< - "npoints="<Next()) { + + Int_t idetector = rawStream->GetDet(); // current detector + Int_t imcm = rawStream->GetMCM(); // current MCM + Int_t irob = rawStream->GetROB(); // current ROB + + //printf("Detector %d\n",idetector); + + if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){ + + // Fill + withInput = TMath::Max(FillDAQ(phvalue),withInput); + + + // reset + for(Int_t k = 0; k < 36; k++){ + for(Int_t j = 0; j < 16; j++){ + for(Int_t c = 0; c < 144; c++){ + phvalue[j][c][k] = 0.0; + } + } + } + } + + fDetectorPreviousTrack = idetector; + fMCMPrevious = imcm; + fROBPrevious = irob; + + nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data + if(nbtimebin == 0) return 0; + if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0; + fTimeMax = nbtimebin; + + //baseline = rawStream->GetCommonAdditive(); // common additive baseline + fNumberClustersf = fTimeMax; + fNumberClusters = (Int_t)(0.6*fTimeMax); + + + Int_t *signal = rawStream->GetSignals(); // current ADC signal + Int_t col = rawStream->GetCol(); + Int_t row = rawStream->GetRow(); + + + //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline); + + + for(Int_t itime = 0; itime < nbtimebin; itime++){ + phvalue[row][col][itime] = signal[itime]-baseline; + } + } - (* fDebugStreamer) << "VDRIFT"<< - "snpright="<= 0.1) return kFALSE; - if(crossrow == 1) return kFALSE; - - if(fLinearFitterOn){ - //Add to the linear fitter of the detector - if( TMath::Abs(snp) < 1.){ - Double_t x = tnp-dzdx*tnt; - (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt); - if(fLinearFitterDebugOn) { - fLinearVdriftFit->Update(detector,x,pars[1]); + // Fill + withInput = TMath::Max(FillDAQ(phvalue),withInput); + + // reset + for(Int_t k = 0; k < 36; k++){ + for(Int_t j = 0; j < 16; j++){ + for(Int_t c = 0; c < 144; c++){ + phvalue[j][c][k] = 0.0; + } + } } - fEntriesLinearFitter[detector]++; } + } - //AliInfo("End of FindP1TrackPH with success!") - return kTRUE; + else{ -} -//____________Offine tracking in the AliTRDtracker_____________________________ -Bool_t AliTRDCalibraFillHisto::HandlePRF() -{ - // - // For the offline tracking - // Fit the tracklet with a line and take the position as reference for the PRF - // + while (rawStream->Next()) { - //Number of points - Int_t npoints = fListClusters->GetEntriesFast(); // number of total points - Int_t nb3pc = 0; // number of three pads clusters used for fit - Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); // detector - + Int_t idetector = rawStream->GetDet(); // current detector + Int_t imcm = rawStream->GetMCM(); // current MCM + Int_t irob = rawStream->GetROB(); // current ROB - // To see the difference due to the fit - Double_t *padPositions; - padPositions = new Double_t[npoints]; - for(Int_t k = 0; k < npoints; k++){ - padPositions[k] = 0.0; - } + //printf("Detector %d\n",idetector); + if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){ - //Find the position by a fit - TLinearFitter fitter(2,"pol1"); - fitter.StoreData(kFALSE); - fitter.ClearPoints(); - for(Int_t k = 0; k < npoints; k++){ - //Take the cluster - AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k); - Short_t *signals = cl->GetSignals(); - Double_t time = cl->GetLocalTimeBin(); - //Calculate x if possible - Float_t xcenter = 0.0; - Bool_t echec = kTRUE; - if((time<=7) || (time>=21)) continue; - // Center 3 balanced: position with the center of the pad - if ((((Float_t) signals[3]) > 0.0) && - (((Float_t) signals[2]) > 0.0) && - (((Float_t) signals[4]) > 0.0)) { - echec = kFALSE; - // Security if the denomiateur is 0 - if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) / - ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) { - xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2])))) - / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) - / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4]))))); + // Fill + withInput = TMath::Max(FillDAQ(phvalue),withInput); + + // reset + for(Int_t k = 0; k < 36; k++){ + for(Int_t j = 0; j < 16; j++){ + for(Int_t c = 0; c < 144; c++){ + phvalue[j][c][k] = 0.0; + } + } + } } - else { - echec = kTRUE; + + fDetectorPreviousTrack = idetector; + fMCMPrevious = imcm; + fROBPrevious = irob; + + //baseline = rawStream->GetCommonAdditive(); // common baseline + + fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data + fNumberClustersf = fTimeMax; + fNumberClusters = (Int_t)(0.6*fTimeMax); + Int_t *signal = rawStream->GetSignals(); // current ADC signal + Int_t col = rawStream->GetCol(); + Int_t row = rawStream->GetRow(); + + + //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline); + + for(Int_t itime = 0; itime < fTimeMax; itime++){ + phvalue[row][col][itime] = signal[itime]-baseline; } } - if(TMath::Abs(xcenter) > 0.5) echec = kTRUE; - if(echec) continue; - //if no echec: calculate with the position of the pad - // Position of the cluster - Double_t padPosition = xcenter + cl->GetPadCol(); - padPositions[k] = padPosition; - nb3pc++; - fitter.AddPoint(&time, padPosition,1); - }//clusters loop + + // fill the last one + if(fDetectorPreviousTrack != -1){ + + // Fill + withInput = TMath::Max(FillDAQ(phvalue),withInput); + + // reset + for(Int_t k = 0; k < 36; k++){ + for(Int_t j = 0; j < 16; j++){ + for(Int_t c = 0; c < 144; c++){ + phvalue[j][c][k] = 0.0; + } + } + } + } + } + + return withInput; + +} +//_____________________________________________________________________ +Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck) +{ + // + // Event processing loop - AliRawReader + // Testbeam 2007 version + // - //printf("nb3pc %d, npoints %d\n",nb3pc,npoints); - if(nb3pc < 3) return kFALSE; - fitter.Eval(); - TVectorD line(2); - fitter.GetParameters(line); - Float_t pointError = -1.0; - pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc); + AliTRDrawStreamBase rawStream(rawReader); + + rawReader->Select("TRD"); + + return ProcessEventDAQ(&rawStream, nocheck); +} + +//_________________________________________________________________________ +Int_t AliTRDCalibraFillHisto::ProcessEventDAQ( +#ifdef ALI_DATE + const eventHeaderStruct *event, + Bool_t nocheck +#else + const eventHeaderStruct* /*event*/, + Bool_t /*nocheck*/ + +#endif + ) +{ + // + // process date event + // Testbeam 2007 version + // +#ifdef ALI_DATE + AliRawReader *rawReader = new AliRawReaderDate((void*)event); + Int_t result=ProcessEventDAQ(rawReader, nocheck); + delete rawReader; + return result; +#else + Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE"); + return 0; +#endif + +} +////////////////////////////////////////////////////////////////////////////// +// Routine inside the DAQ process +///////////////////////////////////////////////////////////////////////////// +//_______________________________________________________________________ +Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){ + + // + // Look for the maximum by collapsing over the time + // Sum over four pad col and two pad row + // + + Int_t used = 0; + + + Int_t idect = fDetectorPreviousTrack; + //printf("Enter Detector %d\n",fDetectorPreviousTrack); + Double_t sum[36]; + for(Int_t tb = 0; tb < 36; tb++){ + sum[tb] = 0.0; + } + + //fGoodTracklet = kTRUE; + //fDetectorPreviousTrack = 0; + + + /////////////////////////// + // look for maximum + ///////////////////////// + + Int_t imaxRow = 0; + Int_t imaxCol = 0; + Double_t integralMax = -1; + for (Int_t ir = 1; ir <= 15; ir++) + { + for (Int_t ic = 2; ic <= 142; ic++) + { + Double_t integral = 0; + for (Int_t ishiftR = 0; ishiftR < 2; ishiftR++) + { + for (Int_t ishiftC = -2; ishiftC < 2; ishiftC++) + { + if (ir + ishiftR >= 1 && ir + ishiftR <= 16 && + ic + ishiftC >= 1 && ic + ishiftC <= 144) + { + + for(Int_t tb = 0; tb< fTimeMax; tb++){ + integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb]; + }// addtb + } //addsignal + } //shiftC + } // shiftR + + if (integralMax < integral) + { + imaxRow = ir; + imaxCol = ic; + integralMax = integral; + } // check max integral + } //ic + } // ir + + //printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax); + + if((imaxRow == 0) || (imaxCol == 0)) { + used=1; + return used; + } + //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol); + //if(!fGoodTracklet) used = 1;; + + // ///////////////////////////////////////////////////// + // sum ober 2 row and 4 pad cols for each time bins + // //////////////////////////////////////////////////// + + + for (Int_t ir = imaxRow - 1; ir < imaxRow + 1; ir++) + { + for (Int_t ic = imaxCol - 2; ic < imaxCol + 2; ic++) + { + for(Int_t it = 0; it < fTimeMax; it++){ + sum[it] += phvalue[ir][ic][it]; + } + }//ic + }//ir - // Now fill the PRF - for(Int_t k = 0; k < npoints; k++){ - //Take the cluster - AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k); - Short_t *signals = cl->GetSignals(); // signal - Double_t time = cl->GetLocalTimeBin(); // time bin - Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit - Float_t padPos = cl->GetPadCol(); // middle pad - Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit - Float_t ycenter = 0.0; // relative center charge - Float_t ymin = 0.0; // relative left charge - Float_t ymax = 0.0; // relative right charge - Double_t tgl = fPar3[(Int_t)time]; // dz/dl and not dz/dx - Double_t pt = fPar4[(Int_t)time]; // pt - Float_t dzdx = 0.0; // dzdx + Int_t nbcl = 0; + Double_t sumcharge = 0.0; + for(Int_t it = 0; it < fTimeMax; it++){ + sumcharge += sum[it]; + if(sum[it] > 20.0) nbcl++; + } - //Requiere simply two pads clusters at least - if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) || - ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){ - Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]); - if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum; - if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum; - if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum; + ///////////////////////////////////////////////////////// + // Debug + //////////////////////////////////////////////////////// + if(fDebugLevel > 0){ + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + Double_t amph0 = sum[0]; + Double_t amphlast = sum[fTimeMax-1]; + Double_t rms = TMath::RMS(fTimeMax,sum); + Int_t goodtracklet = (Int_t) fGoodTracklet; + for(Int_t it = 0; it < fTimeMax; it++){ + Double_t clustera = sum[it]; + + (* fDebugStreamer) << "FillDAQa"<< + "ampTotal="<GetXbins(2)+ grouplocal; // calcul the corresponding group - Double_t snp = fPar2[(Int_t)time]; // sin angle in xy plan - Float_t xcl = cl->GetY(); // y cluster - Float_t qcl = cl->GetQ(); // charge cluster - Int_t plane = GetPlane(detector); // plane - Int_t chamber = GetChamber(detector); // chamber - Double_t xdiff = dpad; // reconstructed position constant - Double_t x = dpad; // reconstructed position moved - Float_t ep = pointError; // error of fit - Float_t signal1 = (Float_t)signals[1]; // signal at the border - Float_t signal3 = (Float_t)signals[3]; // signal - Float_t signal2 = (Float_t)signals[2]; // signal - Float_t signal4 = (Float_t)signals[4]; // signal - Float_t signal5 = (Float_t)signals[5]; // signal at the border - Float_t tnp = 0.0; - if(TMath::Abs(snp) < 1.0){ - tnp = snp / (TMath::Sqrt(1-snp*snp)); - dzdx = tgl*TMath::Sqrt(1+tnp*tnp); + } + + //////////////////////////////////////////////////////// + // fill + /////////////////////////////////////////////////////// + if(sum[0] > 100.0) used = 1; + if(nbcl < fNumberClusters) used = 1; + if(nbcl > fNumberClustersf) used = 1; + + //if(fDetectorPreviousTrack == 15){ + // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]); + //} + //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1; + if(used == 0){ + for(Int_t it = 0; it < fTimeMax; it++){ + if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax); + else{ + if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax); + } } + + //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack)); + used = 2; + //printf("Pass Detector %d\n",fDetectorPreviousTrack); - if(fDebugLevel > 0){ - if ( !fDebugStreamer ) { - //debug stream - TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root"); - if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer - } - - (* fDebugStreamer) << "PRF0"<< - "caligroup="<Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal); + + + return kTRUE; + +} +//____________Write_____________________________________________________ +//_____________________________________________________________________ +void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append) +{ + // + // Write infos to file + // + + //For debugging + if ( fDebugStreamer ) { + delete fDebugStreamer; + fDebugStreamer = 0x0; + } + + AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d" + ,fNumberTrack + ,fNumberUsedCh[0] + ,fNumberUsedCh[1] + ,fNumberUsedPh[0] + ,fNumberUsedPh[1])); + + TDirectory *backup = gDirectory; + TString option; + + if ( append ) + option = "update"; + else + option = "recreate"; + + TFile f(filename,option.Data()); + + TStopwatch stopwatch; + stopwatch.Start(); + if(fVector2d) { + f.WriteTObject(fCalibraVector); + } + + if (fCH2dOn ) { + if (fHisto2d) { + f.WriteTObject(fCH2d); + } + } + if (fPH2dOn ) { + if (fHisto2d) { + f.WriteTObject(fPH2d); + } + } + if (fPRF2dOn) { + if (fHisto2d) { + f.WriteTObject(fPRF2d); } + } + if(fLinearFitterOn){ + AnalyseLinearFitter(); + f.WriteTObject(fLinearVdriftFit); + } + + f.Close(); + + if ( backup ) backup->cd(); + + AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs" + ,stopwatch.RealTime(),stopwatch.CpuTime())); +} +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Stats stuff +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +//___________________________________________probe the histos__________________________________________________ +Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i) +{ + // + // Check the number of stats in h, 0 is TH2I 1 is TProfile2D + // debug mode with 2 for TH2I and 3 for TProfile2D + // It gives a pointer to a Double_t[7] with the info following... + // [0] : number of calibration groups with entries + // [1] : minimal number of entries found + // [2] : calibration group number of the min + // [3] : maximal number of entries found + // [4] : calibration group number of the max + // [5] : mean number of entries found + // [6] : mean relative error + // + + Double_t *info = new Double_t[7]; + + // Number of Xbins (detectors or groups of pads) + Int_t nbins = h->GetNbinsY(); //number of calibration groups + Int_t nxbins = h->GetNbinsX(); //number of bins per histo + + // Initialise + Double_t nbwe = 0; //number of calibration groups with entries + Double_t minentries = 0; //minimal number of entries found + Double_t maxentries = 0; //maximal number of entries found + Double_t placemin = 0; //calibration group number of the min + Double_t placemax = -1; //calibration group number of the max + Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry + Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D + + Double_t counter = 0; + + //Debug + TH1F *nbEntries = 0x0;//distribution of the number of entries + TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group + TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule - // some cuts - if(npoints < fNumberClusters) continue; - if(nb3pc <= 5) continue; - if((time >= 21) || (time < 7)) continue; - if(TMath::Abs(snp) >= 1.0) continue; - if(qcl < 80) continue; + // Beginning of the loop over the calibration groups + for (Int_t idect = 0; idect < nbins; idect++) { + + TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e"); + projch->SetDirectory(0); - Bool_t echec = kFALSE; - Double_t shift = 0.0; - //Calculate the shift in x coresponding to this tnp - if(fNgroupprf != 0.0){ - shift = -3.0*(fNgroupprf-1)-1.5; - Double_t limithigh = -0.2*(fNgroupprf-1); - if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE; - else{ - while(tnp > limithigh){ - limithigh += 0.2; - shift += 3.0; - } + // Number of entries for this calibration group + Double_t nentries = 0.0; + if((i%2) == 0){ + for (Int_t k = 0; k < nxbins; k++) { + nentries += h->GetBinContent(h->GetBin(k+1,idect+1)); } } - if (fHisto2d && !echec) { - if(TMath::Abs(dpad) < 1.5) { - fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter); - fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter); - } - if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) { - fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin); - fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin); - } - if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) { - fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax); - fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax); + else{ + for (Int_t k = 0; k < nxbins; k++) { + nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1)); + if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) { + meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1))))); + counter++; + } } } - //Not equivalent anymore here! - if (fVector2d && !echec) { - if(TMath::Abs(dpad) < 1.5) { - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter); - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter); - } - if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) { - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin); - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin); + + //Debug + if(i > 1){ + if((!((Bool_t)nbEntries)) && (nentries > 0)){ + nbEntries = new TH1F("Number of entries","Number of entries" + ,100,(Int_t)nentries/2,nentries*2); + nbEntries->SetDirectory(0); + nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group" + ,nbins,0,nbins); + nbEntriesPerGroup->SetDirectory(0); + nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule" + ,(Int_t)(nbins/18),0,(Int_t)(nbins/18)); + nbEntriesPerSp->SetDirectory(0); } - if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) { - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax); - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax); + if(nbEntries){ + if(nentries > 0) nbEntries->Fill(nentries); + nbEntriesPerGroup->Fill(idect+0.5,nentries); + nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries); } } + + //min amd max + if(nentries > maxentries){ + maxentries = nentries; + placemax = idect; + } + if(idect == 0) { + minentries = nentries; + } + if(nentries < minentries){ + minentries = nentries; + placemin = idect; + } + //nbwe + if(nentries > 0) { + nbwe++; + meanstats += nentries; + } + }//calibration groups loop + + if(nbwe > 0) meanstats /= nbwe; + if(counter > 0) meanrelativerror /= counter; + + AliInfo(Form("There are %f calibration groups with entries",nbwe)); + AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin)); + AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax)); + AliInfo(Form("The mean number of entries is %f",meanstats)); + if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror)); + + info[0] = nbwe; + info[1] = minentries; + info[2] = placemin; + info[3] = maxentries; + info[4] = placemax; + info[5] = meanstats; + info[6] = meanrelativerror; + + if(i > 1){ + gStyle->SetPalette(1); + gStyle->SetOptStat(1111); + gStyle->SetPadBorderMode(0); + gStyle->SetCanvasColor(10); + gStyle->SetPadLeftMargin(0.13); + gStyle->SetPadRightMargin(0.01); + TCanvas *stat = new TCanvas("stat","",50,50,600,800); + stat->Divide(2,1); + stat->cd(1); + nbEntries->Draw(""); + stat->cd(2); + nbEntriesPerSp->SetStats(0); + nbEntriesPerSp->Draw(""); + TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800); + stat1->cd(); + nbEntriesPerGroup->SetStats(0); + nbEntriesPerGroup->Draw(""); } - return kTRUE; - + + return info; + } -//____________Offine tracking in the AliTRDtracker_____________________________ -Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrack(AliTRDtrack *t, Int_t index0, Int_t index1) +//____________________________________________________________________________ +Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH() { // - // For the offline tracking - // This function will be called in the functions UpdateHistogram... - // to fill the find the parameter P1 of a track for the drift velocity calibration + // Return a Int_t[4] with: + // 0 Mean number of entries + // 1 median of number of entries + // 2 rms of number of entries + // 3 number of group with entries // - - //Number of points: if less than 3 return kFALSE - Int_t npoints = index1-index0; - if(npoints <= 2) return kFALSE; - - //Variables - TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet - Double_t snp = 0.0; // sin angle in the plan yx track - Double_t y = 0.0; // y clusters in the middle of the chamber - Double_t z = 0.0; // z cluster in the middle of the chamber - Double_t dydt = 0.0; // dydt tracklet after straight line fit - Double_t tnp = 0.0; // tan angle in the plan xy track - Double_t tgl = 0.0; // dz/dl and not dz/dx! - Double_t errorpar = 0.0; // error after straight line fit on dy/dt - Double_t pointError = 0.0; // error after straight line fit - Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector - //Int_t snpright = 1; // if we took in the middle snp - Int_t crossrow = 0; // if it crosses a pad row - Double_t tiltingangle = 0; // tiltingangle of the pad - Float_t dzdx = 0; // dz/dx now from dz/dl - Int_t nbli = 0; // number linear fitter points - AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector)); - - linearFitterTracklet.StoreData(kFALSE); - linearFitterTracklet.ClearPoints(); - - //if more than one row - Int_t rowp = -1; // if it crosses a pad row - - //tiltingangle - tiltingangle = padplane->GetTiltingAngle(); - Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle + Double_t *stat = new Double_t[4]; + stat[3] = 0.0; - //Fill with points - for(Int_t k = 0; k < npoints; k++){ - - AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0); - Double_t ycluster = cl->GetY(); - Int_t time = cl->GetLocalTimeBin(); - Double_t timeis = time/fSf; - //See if cross two pad rows - Int_t row = padplane->GetPadRowNumber(cl->GetZ()); - if(k==0) rowp = row; - if(row != rowp) crossrow = 1; - //Take in the middle of the chamber - //FollowBack - //if(time > (Int_t) 10) { - //Follow - if(time < (Int_t) 11) { - z = cl->GetZ(); - y = cl->GetY(); + Int_t nbofgroups = CalculateTotalNumberOfBins(0); + Double_t *weight = new Double_t[nbofgroups]; + Int_t *nonul = new Int_t[nbofgroups]; + + for(Int_t k = 0; k < nbofgroups; k++){ + if(fEntriesCH[k] > 0) { + weight[k] = 1.0; + nonul[(Int_t)stat[3]] = fEntriesCH[k]; + stat[3]++; } - linearFitterTracklet.AddPoint(&timeis,ycluster,1); - nbli++; + else weight[k] = 0.0; } + stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight); + stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight); + stat[2] = TMath::RMS((Int_t)stat[3],nonul); - // take now the snp, tnp and tgl from the track - snp = t->GetSnpPlane(GetPlane(detector)); - tgl = t->GetTglPlane(GetPlane(detector)); + return stat; - //FollowBack - //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() < 10) snpright = 0; - //Follow - //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() >= 11) snpright = 0; - if(nbli <= 2) return kFALSE; - - // Do the straight line fit now - TVectorD pars; - linearFitterTracklet.Eval(); - linearFitterTracklet.GetParameters(pars); - pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli); - errorpar = linearFitterTracklet.GetParError(1)*pointError; - dydt = pars[1]; - - if( TMath::Abs(snp) < 1.){ - tnp = snp / (TMath::Sqrt(1-(snp*snp))); - } - dzdx = tgl*TMath::Sqrt(1+tnp*tnp); +} +//____________________________________________________________________________ +Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const +{ + // + // Return a Int_t[4] with: + // 0 Mean number of entries + // 1 median of number of entries + // 2 rms of number of entries + // 3 number of group with entries + // - if(fDebugLevel > 0){ - if ( !fDebugStreamer ) { - //debug stream - TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root"); - if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer - } - - - (* fDebugStreamer) << "VDRIFT"<< - //"snpright="<= 0.1) return kFALSE; - if(crossrow == 1) return kFALSE; - - if(fLinearFitterOn){ - //Add to the linear fitter of the detector - if( TMath::Abs(snp) < 1.){ - Double_t x = tnp-dzdx*tnt; - (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt); - if(fLinearFitterDebugOn) { - fLinearVdriftFit->Update(detector,x,pars[1]); - } - fEntriesLinearFitter[detector]++; + Int_t nbofgroups = 540; + Double_t *weight = new Double_t[nbofgroups]; + Int_t *nonul = new Int_t[nbofgroups]; + + for(Int_t k = 0; k < nbofgroups; k++){ + if(fEntriesLinearFitter[k] > 0) { + weight[k] = 1.0; + nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k]; + stat[3]++; } + else weight[k] = 0.0; } - //AliInfo("End of FindP1TrackPH with success!") - return kTRUE; + stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight); + stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight); + stat[2] = TMath::RMS((Int_t)stat[3],nonul); + + return stat; } -//____________Offine tracking in the AliTRDtracker_____________________________ -Bool_t AliTRDCalibraFillHisto::HandlePRFtrack(AliTRDtrack *t, Int_t index0, Int_t index1) +////////////////////////////////////////////////////////////////////////////////////// +// Create Histos +////////////////////////////////////////////////////////////////////////////////////// +//_____________________________________________________________________________ +void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn) { // - // For the offline tracking - // Fit the tracklet with a line and take the position as reference for the PRF + // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each + // If fNgroupprf is zero then no binning in tnp // - //Number of points - Int_t npoints = index1-index0; // number of total points - Int_t nb3pc = 0; // number of three pads clusters used for fit - Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector - - - // To see the difference due to the fit - Double_t *padPositions; - padPositions = new Double_t[npoints]; - for(Int_t k = 0; k < npoints; k++){ - padPositions[k] = 0.0; - } - - - //Find the position by a fit - TLinearFitter fitter(2,"pol1"); - fitter.StoreData(kFALSE); - fitter.ClearPoints(); - for(Int_t k = 0; k < npoints; k++){ - //Take the cluster - AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0); - Short_t *signals = cl->GetSignals(); - Double_t time = cl->GetLocalTimeBin(); - //Calculate x if possible - Float_t xcenter = 0.0; - Bool_t echec = kTRUE; - if((time<=7) || (time>=21)) continue; - // Center 3 balanced: position with the center of the pad - if ((((Float_t) signals[3]) > 0.0) && - (((Float_t) signals[2]) > 0.0) && - (((Float_t) signals[4]) > 0.0)) { - echec = kFALSE; - // Security if the denomiateur is 0 - if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) / - ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) { - xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2])))) - / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) - / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4]))))); - } - else { - echec = kTRUE; - } - } - if(TMath::Abs(xcenter) > 0.5) echec = kTRUE; - if(echec) continue; - //if no echec: calculate with the position of the pad - // Position of the cluster - Double_t padPosition = xcenter + cl->GetPadCol(); - padPositions[k] = padPosition; - nb3pc++; - fitter.AddPoint(&time, padPosition,1); - }//clusters loop + TString name("Nz"); + name += fCalibraMode->GetNz(2); + name += "Nrphi"; + name += fCalibraMode->GetNrphi(2); + name += "Ngp"; + name += fNgroupprf; - //printf("nb3pc %d, npoints %d\n",nb3pc,npoints); - if(nb3pc < 3) return kFALSE; - fitter.Eval(); - TVectorD line(2); - fitter.GetParameters(line); - Float_t pointError = -1.0; - pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc); - - // Take the tgl and snp with the track t now - Double_t tgl = t->GetTglPlane(GetPlane(detector)); //dz/dl and not dz/dx - Double_t snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in xy plan - Float_t dzdx = 0.0; // dzdx - Float_t tnp = 0.0; - if(TMath::Abs(snp) < 1.0){ - tnp = snp / (TMath::Sqrt(1-snp*snp)); - dzdx = tgl*TMath::Sqrt(1+tnp*tnp); + if(fNgroupprf != 0){ + + fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name + ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn); + fPRF2d->SetYTitle("Det/pad groups"); + fPRF2d->SetXTitle("Position x/W [pad width units]"); + fPRF2d->SetZTitle("Q_{i}/Q_{total}"); + fPRF2d->SetStats(0); + } + else{ + fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name + ,fNumberBinPRF,-1.5,1.5,nn,0,nn); + fPRF2d->SetYTitle("Det/pad groups"); + fPRF2d->SetXTitle("Position x/W [pad width units]"); + fPRF2d->SetZTitle("Q_{i}/Q_{total}"); + fPRF2d->SetStats(0); } +} - // Now fill the PRF - for(Int_t k = 0; k < npoints; k++){ - //Take the cluster - AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0); - Short_t *signals = cl->GetSignals(); // signal - Double_t time = cl->GetLocalTimeBin(); // time bin - Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit - Float_t padPos = cl->GetPadCol(); // middle pad - Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit - Float_t ycenter = 0.0; // relative center charge - Float_t ymin = 0.0; // relative left charge - Float_t ymax = 0.0; // relative right charge +//_____________________________________________________________________________ +void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn) +{ + // + // Create the 2D histos + // + + TString name("Nz"); + name += fCalibraMode->GetNz(1); + name += "Nrphi"; + name += fCalibraMode->GetNrphi(1); + fPH2d = new TProfile2D("PH2d",(const Char_t *) name + ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf + ,nn,0,nn); + fPH2d->SetYTitle("Det/pad groups"); + fPH2d->SetXTitle("time [#mus]"); + fPH2d->SetZTitle(" [a.u.]"); + fPH2d->SetStats(0); +} +//_____________________________________________________________________________ +void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn) +{ + // + // Create the 2D histos + // - //Requiere simply two pads clusters at least - if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) || - ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){ - Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]); - if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum; - if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum; - if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum; - } - - //calibration group - Int_t *rowcol = CalculateRowCol(cl); // calcul col and row pad of the cluster - Int_t grouplocal = CalculateCalibrationGroup(2,rowcol); // calcul the corresponding group - Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group - Float_t xcl = cl->GetY(); // y cluster - Float_t qcl = cl->GetQ(); // charge cluster - Int_t plane = GetPlane(detector); // plane - Int_t chamber = GetChamber(detector); // chamber - Double_t xdiff = dpad; // reconstructed position constant - Double_t x = dpad; // reconstructed position moved - Float_t ep = pointError; // error of fit - Float_t signal1 = (Float_t)signals[1]; // signal at the border - Float_t signal3 = (Float_t)signals[3]; // signal - Float_t signal2 = (Float_t)signals[2]; // signal - Float_t signal4 = (Float_t)signals[4]; // signal - Float_t signal5 = (Float_t)signals[5]; // signal at the border - - + TString name("Nz"); + name += fCalibraMode->GetNz(0); + name += "Nrphi"; + name += fCalibraMode->GetNrphi(0); - if(fDebugLevel > 0){ - if ( !fDebugStreamer ) { - //debug stream - TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root"); - if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer - } - - (* fDebugStreamer) << "PRF0"<< - "caligroup="<= 21) || (time < 7)) continue; - if(TMath::Abs(snp) >= 1.0) continue; - if(qcl < 80) continue; - - Bool_t echec = kFALSE; - Double_t shift = 0.0; - //Calculate the shift in x coresponding to this tnp - if(fNgroupprf != 0.0){ - shift = -3.0*(fNgroupprf-1)-1.5; - Double_t limithigh = -0.2*(fNgroupprf-1); - if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE; - else{ - while(tnp > limithigh){ - limithigh += 0.2; - shift += 3.0; - } - } - } - if (fHisto2d && !echec) { - if(TMath::Abs(dpad) < 1.5) { - fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter); - fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter); - } - if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) { - fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin); - fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin); - } - if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) { - fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax); - fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax); - } - } - //Not equivalent anymore here! - if (fVector2d && !echec) { - if(TMath::Abs(dpad) < 1.5) { - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter); - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter); - } - if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) { - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin); - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin); - } - if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) { - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax); - fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax); - } - } + fCH2d = new TH2I("CH2d",(const Char_t *) name + ,fNumberBinCharge,0,300,nn,0,nn); + fCH2d->SetYTitle("Det/pad groups"); + fCH2d->SetXTitle("charge deposit [a.u]"); + fCH2d->SetZTitle("counts"); + fCH2d->SetStats(0); + fCH2d->Sumw2(); + +} +////////////////////////////////////////////////////////////////////////////////// +// Set relative scale +///////////////////////////////////////////////////////////////////////////////// +//_____________________________________________________________________________ +void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale) +{ + // + // Set the factor that will divide the deposited charge + // to fit in the histo range [0,300] + // + + if (RelativeScale > 0.0) { + fRelativeScale = RelativeScale; + } + else { + AliInfo("RelativeScale must be strict positif!"); } - return kTRUE; + +} +////////////////////////////////////////////////////////////////////////////////// +// Quick way to fill a histo +////////////////////////////////////////////////////////////////////////////////// +//_____________________________________________________________________ +void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y) +{ + // + // FillCH2d: Marian style + // + + //skip simply the value out of range + if((y>=300.0) || (y<0.0)) return; + + //Calcul the y place + Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1; + Int_t place = (fNumberBinCharge+2)*(x+1)+yplace; + //Fill + fCH2d->GetArray()[place]++; + } -// -//____________Some basic geometry function_____________________________________ -// + +////////////////////////////////////////////////////////////////////////////////// +// Geometrical functions +/////////////////////////////////////////////////////////////////////////////////// //_____________________________________________________________________________ -Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const +Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const { // - // Reconstruct the plane number from the detector number + // Reconstruct the layer number from the detector number // return ((Int_t) (d % 6)); @@ -2959,14 +3080,14 @@ Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const } //_____________________________________________________________________________ -Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const +Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const { // - // Reconstruct the chamber number from the detector number + // Reconstruct the stack number from the detector number // - Int_t fgkNplan = 6; + const Int_t kNlayer = 6; - return ((Int_t) (d % 30) / fgkNplan); + return ((Int_t) (d % 30) / kNlayer); } @@ -2981,6 +3102,9 @@ Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const return ((Int_t) (d / fg)); } +/////////////////////////////////////////////////////////////////////////////////// +// Getter functions for DAQ of the CH2d and the PH2d +////////////////////////////////////////////////////////////////////////////////// //_____________________________________________________________________ TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency) { @@ -2995,31 +3119,28 @@ TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequ fTimeMax = nbtimebin; fSf = samplefrequency; - /* - AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d" - ,fCalibraMode->GetNz(1) - ,fCalibraMode->GetNrphi(1))); - - // Calcul the number of Xbins - Int_t ntotal1 = 0; - fCalibraMode->ModePadCalibration(2,1); - fCalibraMode->ModePadFragmentation(0,2,0,1); - fCalibraMode->SetDetChamb2(1); - ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1); - fCalibraMode->ModePadCalibration(0,1); - fCalibraMode->ModePadFragmentation(0,0,0,1); - fCalibraMode->SetDetChamb0(1); - ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1); - AliInfo(Form("Total number of Xbins: %d",ntotal1)); - - CreatePH2d(ntotal1); - */ - CreatePH2d(540); return fPH2d; } //_____________________________________________________________________ +TH2I* AliTRDCalibraFillHisto::GetCH2d() +{ + // + // return pointer to fCH2d TH2I + // create a new TH2I if it doesn't exist allready + // + if ( fCH2d ) + return fCH2d; + + CreateCH2d(540); + + return fCH2d; +} +//////////////////////////////////////////////////////////////////////////////////////////// +// Drift velocity calibration +/////////////////////////////////////////////////////////////////////////////////////////// +//_____________________________________________________________________ TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force) { // @@ -3039,25 +3160,6 @@ TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t fo return linearfitter; } -//_____________________________________________________________________ -void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y) -{ - // - // FillCH2d: Marian style - // - - //skip simply the value out of range - if((y>=300.0) || (y<0.0)) return; - - //Calcul the y place - Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1; - Int_t place = (fNumberBinCharge+2)*(x+1)+yplace; - - //Fill - fCH2d->GetArray()[place]++; - -} - //____________________________________________________________________________ void AliTRDCalibraFillHisto::AnalyseLinearFitter() { @@ -3084,20 +3186,3 @@ void AliTRDCalibraFillHisto::AnalyseLinearFitter() } } } -//________________________________________________________________________________ -Int_t AliTRDCalibraFillHisto::Arrondi(Double_t x) const -{ - // Partie entiere of the (x+0.5) - - int i; - if (x >= (-0.5)) { - i = int(x + 0.5); - //if (x + 0.5 == Float_t(i) && i & 1) i--; - } else { - i = int(x - 0.5); - //if (x - 0.5 == Float_t(i) && i & 1) i++; - if((x-0.5)==i) i++; - - } - return i; -}