From: cblume Date: Fri, 11 May 2007 17:40:29 +0000 (+0000) Subject: Go back to previous PID method X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=b08a2e7e25bbcc40be9264dcce71c6635cff8a80;p=u%2Fmrichter%2FAliRoot.git Go back to previous PID method --- diff --git a/TRD/AliTRDtracker.cxx b/TRD/AliTRDtracker.cxx index 9efdabc4de9..a180ce49322 100644 --- a/TRD/AliTRDtracker.cxx +++ b/TRD/AliTRDtracker.cxx @@ -58,9 +58,9 @@ ClassImp(AliTRDtracker) -const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5; -const Float_t AliTRDtracker::fgkLabelFraction = 0.8; // ?? -const Double_t AliTRDtracker::fgkMaxChi2 = 12.0; +const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5; // +const Float_t AliTRDtracker::fgkLabelFraction = 0.8; // +const Double_t AliTRDtracker::fgkMaxChi2 = 12.0; // const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // Corresponds to tan = 3 const Double_t AliTRDtracker::fgkMaxStep = 2.0; // Maximal step size in propagation @@ -3468,65 +3468,57 @@ Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c) void AliTRDtracker::CookdEdxTimBin(AliTRDtrack &TRDtrack) { // - // Build the information needed for TRD PID calculations. This are: - // - dE/dx - total - // - dE/dx - slices (0-14, 15-22, NONE). The Edep per slice is - // calculated as the sum of clusters charge weighted with average - // value from the spectrum - // - Time bin for Max. Cluster - // - // Authors: + // This is setting fdEdxPlane and fTimBinPlane + // Sums up the charge in each plane for track TRDtrack and also get the + // Time bin for Max. Cluster // Prashant Shukla (shukla@physi.uni-heidelberg.de) - // Alex Bercuci (a.bercuci@gsi.de) + // Double_t clscharge[AliESDtrack::kNPlane][AliESDtrack::kNSlice]; Double_t maxclscharge[AliESDtrack::kNPlane]; Int_t nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice]; Int_t timebin[AliESDtrack::kNPlane]; - // Initialization of variables + // Initialization of cluster charge per plane. for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) { for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) { - clscharge[iPlane][iSlice] = 0.; + clscharge[iPlane][iSlice] = 0.0; nCluster[iPlane][iSlice] = 0; } + } + + // Initialization of cluster charge per plane. + for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) { timebin[iPlane] = -1; maxclscharge[iPlane] = 0.0; } - // Average PH spectrum used for weighting the edep. It has to - // come eather from caligration DB or as parameter of the class. To - // be discuss with software responsibles - const Float_t fAveragePH[] = { 5.010e-03, 2.017e-02, 6.221e-02, 6.706e-02 - , 5.551e-02, 4.812e-02, 4.395e-02, 4.219e-02 - , 4.172e-02, 4.156e-02, 4.162e-02, 4.204e-02 - , 4.278e-02, 4.367e-02, 4.460e-02, 4.554e-02 - , 4.674e-02, 4.810e-02, 5.005e-02, 5.258e-02 - , 5.587e-02, 5.892e-02, 5.587e-02, 5.258e-02 }; - // Loop through all clusters associated to track TRDtrack - AliTRDcluster *pTRDcluster = 0x0; - for (Int_t iClus = 0; iClus < TRDtrack.GetNumberOfClusters(); iClus++) { + Int_t nClus = TRDtrack.GetNumberOfClusters(); // from Kalmantrack + for (Int_t iClus = 0; iClus < nClus; iClus++) { Double_t charge = TRDtrack.GetClusterdQdl(iClus); Int_t index = TRDtrack.GetClusterIndex(iClus); - if (!(pTRDcluster = (AliTRDcluster*)GetCluster(index))) continue; - - // check negative time bins ?!?!? - Int_t tb = pTRDcluster->GetLocalTimeBin(); - - // temporary until fix in GetT0 - if(tb<0){ - Warning("CookdEdxTimBin()", Form("Negative time bin value in track %d cluster %d", TRDtrack.GetLabel(), iClus)); - continue; - } - + AliTRDcluster *pTRDcluster = (AliTRDcluster *) GetCluster(index); + if (!pTRDcluster) { + continue; + } + Int_t tb = pTRDcluster->GetLocalTimeBin(); + if (!tb) { + continue; + } Int_t detector = pTRDcluster->GetDetector(); Int_t iPlane = fGeom->GetPlane(detector); - - // 2 slices method - Int_t iSlice = (tb < AliTRDtrack::kMidTimeBin) ? 0 : 1; - // weighted edep with - clscharge[iPlane][iSlice] += charge * fAveragePH[tb]; + if (iPlane >= AliESDtrack::kNPlane) { + AliError(Form("Wrong plane %d",iPlane)); + continue; + } + Int_t iSlice = tb * AliESDtrack::kNSlice + / AliTRDcalibDB::Instance()->GetNumberOfTimeBins(); + if (iSlice >= AliESDtrack::kNSlice) { + AliError(Form("Wrong slice %d",iSlice)); + continue; + } + clscharge[iPlane][iSlice] = clscharge[iPlane][iSlice] + charge; if (charge > maxclscharge[iPlane]) { maxclscharge[iPlane] = charge; timebin[iPlane] = tb; @@ -3534,18 +3526,31 @@ void AliTRDtracker::CookdEdxTimBin(AliTRDtrack &TRDtrack) nCluster[iPlane][iSlice]++; } // End of loop over cluster - // Setting the fdEdxPlane and fTimBinPlane variables - Double_t totalCharge = 0.; + // Setting the fdEdxPlane and fTimBinPlane variabales + Double_t totalCharge = 0.0; + for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) { - for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) { - //if (nCluster[iPlane][iSlice]) { - // clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice]; - //} + for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) { + if (nCluster[iPlane][iSlice]) { + clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice]; + } TRDtrack.SetPIDsignals(clscharge[iPlane][iSlice],iPlane,iSlice); - totalCharge += clscharge[iPlane][iSlice]; + totalCharge = totalCharge+clscharge[iPlane][iSlice]; } - TRDtrack.SetPIDTimBin(timebin[iPlane],iPlane); + TRDtrack.SetPIDTimBin(timebin[iPlane],iPlane); } + + // Still needed ???? + // Int_t i; + // Int_t nc=TRDtrack.GetNumberOfClusters(); + // Float_t dedx=0; + // for (i=0; iGetXaxis(); nbinsx = ax->GetNbins(); ay = hist->GetYaxis(); nbinsy = ay->GetNbins(); - Bool_t kX = (dedx[0] < ax->GetBinUpEdge(nbinsx)); - Bool_t kY = (dedx[1] < ay->GetBinUpEdge(nbinsy)); + Float_t x = dedx[0]+dedx[1], y = dedx[2]; + Bool_t kX = (x < ax->GetBinUpEdge(nbinsx)); + Bool_t kY = (y < ay->GetBinUpEdge(nbinsy)); if(kX) - if(kY) LQ1 = hist->GetBinContent( hist->FindBin(dedx[0], dedx[1])); //fEstimator->Estimate2D2(hist, (Float_t&)dedx[0], (Float_t&)dedx[1]); - else LQ1 = hist->GetBinContent(ax->FindBin(dedx[0]), nbinsy); + if(kY) LQ1 = hist->GetBinContent( hist->FindBin(x, y)); + //fEstimator->Estimate2D2(hist, x, y); + else LQ1 = hist->GetBinContent(ax->FindBin(x), nbinsy); else - if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(dedx[1])); + if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(y)); else LQ1 = hist->GetBinContent(nbinsx, nbinsy); @@ -436,10 +438,11 @@ Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Double_t mom, Double_t *dedx return LQ1; } if(kX) - if(kY) LQ2 = hist->GetBinContent( hist->FindBin(dedx[0], dedx[1])); //fEstimator->Estimate2D2(hist, (Float_t&)dedx[0], (Float_t&)dedx[1]); - else LQ2 = hist->GetBinContent(ax->FindBin(dedx[0]), nbinsy); + if(kY) LQ2 = hist->GetBinContent( hist->FindBin(x, y)); + //fEstimator->Estimate2D2(hist, x, y); + else LQ2 = hist->GetBinContent(ax->FindBin(x), nbinsy); else - if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(dedx[1])); + if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(y)); else LQ2 = hist->GetBinContent(nbinsx, nbinsy);