o bug fixes for 11h gain scenarios (Mikolaj)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Nov 2012 16:07:24 +0000 (16:07 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Nov 2012 16:07:24 +0000 (16:07 +0000)
o better treatment of eta maps for the 'tuned on data' setting (Benjamin)
o load default splines for pPb MC
J. Wiechula

STEER/AOD/AliAODpidUtil.cxx
STEER/ESD/AliESDpid.cxx
STEER/STEERBase/AliPIDResponse.cxx
STEER/STEERBase/AliTPCPIDResponse.cxx
STEER/STEERBase/AliTPCPIDResponse.h

index d272d78..b758b09 100644 (file)
@@ -44,8 +44,6 @@ Float_t AliAODpidUtil::GetTPCsignalTunedOnData(const AliVTrack *t) const {
     Float_t dedx = track->GetTPCsignalTunedOnData();
     if(dedx > 0) return dedx;
 
-    Double_t mom = t->GetTPCmomentum();
-
     dedx = t->GetTPCsignal();
     track->SetTPCsignalTunedOnData(dedx);
 
@@ -92,10 +90,11 @@ Float_t AliAODpidUtil::GetTPCsignalTunedOnData(const AliVTrack *t) const {
            kGood = kFALSE;
 
        if(kGood){
-           Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
-           Double_t sigma=fTPCResponse.GetExpectedSigma(mom,t->GetTPCsignalN(),type);
-           dedx = gRandom->Gaus(bethe,sigma);
-
+           //TODO maybe introduce different dEdxSources?
+        Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
+        Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
+        dedx = gRandom->Gaus(bethe,sigma);
+        
            if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
        }
 
index 4c5df6f..19548d8 100644 (file)
@@ -75,7 +75,6 @@ Float_t AliESDpid::GetTPCsignalTunedOnData(const AliVTrack *t) const {
 
     if(dedx > 0) return dedx;
 
-    Double_t mom = t->GetTPCmomentum();
     dedx = t->GetTPCsignal();
     track->SetTPCsignalTunedOnData(dedx);
     if(dedx < 20) return dedx;
@@ -123,8 +122,9 @@ Float_t AliESDpid::GetTPCsignalTunedOnData(const AliVTrack *t) const {
                kGood = kFALSE;
 
            if(kGood){
-               Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
-               Double_t sigma=fTPCResponse.GetExpectedSigma(mom,t->GetTPCsignalN(),type);
+        //TODO maybe introduce different dEdxSources?
+        Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
+        Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
                dedx = gRandom->Gaus(bethe,sigma);
                if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
            }
index 3dcfb1c..6aacb0c 100644 (file)
@@ -718,7 +718,7 @@ void AliPIDResponse::SetRecoInfo()
   if (fRun >= 186636  ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ }
 
   //exception new pp MC productions from 2011
-  if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
+  if ( (fBeamType=="PP" || fBeamType=="PPB") && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
   // exception for 11f1
   if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
 }
@@ -780,7 +780,7 @@ TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refine
       Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
       Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
       
-      /*TODO NOW NOW
+      /*
       linExtrapolation->ClearPoints();
       
       // For interpolation: Just take the corresponding bin from the old histo.
@@ -847,7 +847,7 @@ TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refine
       Double_t interpolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
                                  + linExtrapolation->GetParameter(2) * centerY;
       */
-      Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;//TODO NOW NOW NOW
+      Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;
       hRefined->SetBinContent(binX, binY, interpolatedValue);      
     }
   } 
@@ -865,48 +865,54 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac
   // Load the TPC eta correction maps from the OADB
   //
   
+  if (fUseTPCEtaCorrection == kFALSE) {
+    // Disable eta correction via setting no maps
+    if (!fTPCResponse.SetEtaCorrMap(0x0))
+      AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled");
+    else
+      AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
+    
+    if (!fTPCResponse.SetSigmaParams(0x0, 0))
+      AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma");
+    else  
+      AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
+    
+    return;
+  }
+
   TString dataType = "DATA";
   TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
   
   if (fIsMC)  {
-    dataType = "MC";
+    if (!fTuneMConData) {
+      period=fMCperiodTPC;
+      dataType="MC";
+    }
     fRecoPass = 1;
     
-    if (fMCperiodTPC.IsNull()) {
+    if (!fTuneMConData && fMCperiodTPC.IsNull()) {
       AliFatal("MC detected, but no MC period set -> Not changing eta maps!");
       return;
     }
-  
-    period = fMCperiodTPC;
   }
+
+  Int_t recopass = fRecoPass;
+  if (fTuneMConData)
+    recopass = fRecoPassUser;
   
-  TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), fRecoPass);
+  TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), recopass);
   
-  AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), fRecoPass));
-  if (fUseTPCEtaCorrection == kFALSE) {
-    // Disable eta correction via setting no maps
-    if (!fTPCResponse.SetEtaCorrMap(0x0))
-      AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled"); 
-    else
-      AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
-    
-    if (!fTPCResponse.SetSigmaParams(0x0, 0))
-      AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma"); 
-    else
-      AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
-    
-    return;
-  }
+  AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), recopass));
   
   // Invalidate old maps
   fTPCResponse.SetEtaCorrMap(0x0);
   fTPCResponse.SetSigmaParams(0x0, 0);
   
   // Load the eta correction maps
-  AliOADBContainer etaMapsCont(Form("TPCetaMaps_%s_pass%d", dataType.Data(), fRecoPass)); 
+  AliOADBContainer etaMapsCont(Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass)); 
   
   Int_t statusCont = etaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
-                                              Form("TPCetaMaps_%s_pass%d", dataType.Data(), fRecoPass));
+                                              Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
   if (statusCont) {
     AliError("Failed initializing TPC eta correction maps from OADB -> Disabled eta correction");
   }
@@ -915,8 +921,8 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac
     
     TH2D* etaMap = 0x0;
     
-    if (fIsMC) {
-      TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), fRecoPass);
+    if (fIsMC && !fTuneMConData) {
+      TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
       etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(searchMap.Data()));
       if (!etaMap) {
         // Try default object
@@ -953,10 +959,10 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac
   }
   
   // Load the sigma parametrisation (1/dEdx vs tanTheta_local (~eta))
-  AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), fRecoPass)); 
+  AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass)); 
   
   statusCont = etaSigmaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
-                                             Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), fRecoPass));
+                                             Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
   if (statusCont) {
     AliError("Failed initializing TPC eta sigma maps from OADB -> Using old sigma parametrisation");
   }
@@ -965,8 +971,8 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac
     
     TObjArray* etaSigmaPars = 0x0;
     
-    if (fIsMC) {
-      TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), fRecoPass);
+    if (fIsMC && !fTuneMConData) {
+      TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
       etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(searchMap.Data()));
       if (!etaSigmaPars) {
         // Try default object
@@ -1096,21 +1102,22 @@ void AliPIDResponse::SetTPCParametrisation()
   TString period=fLHCperiod;
   if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
 
-  AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+  Int_t recopass = fRecoPass;
+  if(fTuneMConData) recopass = fRecoPassUser;
+    
+  AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
   Bool_t found=kFALSE;
   //
   //set the new PID splines
   //
   if (fArrPidResponseMaster){
-    Int_t recopass = fRecoPass;
-    if(fTuneMConData) recopass = fRecoPassUser;
     //for MC don't use period information
     //if (fIsMC) period="[A-Z0-9]*";
     //for MC use MC period information
     //pattern for the default entry (valid for all particles)
     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
 
-    //find particle id ang gain scenario
+    //find particle id and gain scenario
     for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
     {
       TObject *grAll=NULL;
@@ -1177,7 +1184,7 @@ void AliPIDResponse::SetTPCParametrisation()
   else AliInfo("no fArrPidResponseMaster");
 
   if (!found){
-    AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+    AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
   }
 
   //
@@ -1197,7 +1204,7 @@ void AliPIDResponse::SetTPCParametrisation()
   }
   
   if (fArrPidResponseMaster)
-  fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+  fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
   
   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
 
@@ -1736,18 +1743,10 @@ Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID:
   
   Double_t nSigma = -999.;
   
-  //TODO: TPCsignalTunedOnData not taken into account for the new TPCPIDresponse functions, so take the old functions
-  if (fTuneMConData) {
-    Double_t mom  = track->GetTPCmomentum();
-    Double_t sig  = this->GetTPCsignalTunedOnData(track);
-    UInt_t   sigN = track->GetTPCsignalN();
+  if (fTuneMConData)
+    this->GetTPCsignalTunedOnData(track);
   
-    if (sigN>0)
-      nSigma = fTPCResponse.GetNumberOfSigmas(mom, sig, sigN, type);
-  }
-  else {
-    nSigma = fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
-  }
+  nSigma = fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
   
   return nSigma;
 }
@@ -1910,16 +1909,9 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const A
   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
     AliPID::EParticleType type=AliPID::EParticleType(j);
     
-    //TODO: TPCsignalTunedOnData not taken into account for the new TPCPIDresponse functions, so take the old functions
-    if (fTuneMConData) {
-      Double_t mom = track->GetTPCmomentum();
-      bethe=fTPCResponse.GetExpectedSignal(mom,type);
-      sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
-    }
-    else {
-      bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
-      sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
-    }
+    bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+    sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+    
     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
     } else {
index c5daeb5..9b8c515 100644 (file)
@@ -188,12 +188,14 @@ AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
 
   delete fhEtaCorr;
+  fhEtaCorr=0x0;
   if (that.fhEtaCorr){
     fhEtaCorr = new TH2D(*(that.fhEtaCorr));
     fhEtaCorr->SetDirectory(0);
   }
   
   delete fhEtaSigmaPar1;
+  fhEtaSigmaPar1=0x0;
   if (that.fhEtaSigmaPar1){
     fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
     fhEtaSigmaPar1->SetDirectory(0);
@@ -534,7 +536,10 @@ Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
 Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track, 
                                                  AliPID::EParticleType species,
                                                  ETPCdEdxSource dedxSource,
-                                                 Double_t& dEdx, Int_t& nPoints, ETPCgainScenario& gainScenario, TSpline3** responseFunction) const 
+                                                 Double_t& dEdx, 
+                                                 Int_t& nPoints, 
+                                                 ETPCgainScenario& gainScenario, 
+                                                 TSpline3** responseFunction) const 
 {
   // Calculates the right parameters for PID
   //   dEdx parametrization for the proper gain scenario, dEdx 
@@ -547,7 +552,14 @@ Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
   if (dedxSource == kdEdxDefault) {
     // Fast handling for default case. In addition: Keep it simple (don't call additional functions) to
     // avoid possible bugs
-    dEdx = track->GetTPCsignal();
+    
+    // GetTPCsignalTunedOnData will be non-positive, if it has not been set (i.e. in case of MC NOT tuned to data).
+    // If this is the case, just take the normal signal
+    dEdx = track->GetTPCsignalTunedOnData();
+    if (dEdx <= 0) {
+      dEdx = track->GetTPCsignal();
+    }
+    
     nPoints = track->GetTPCsignalN();
     gainScenario = kDefault;
     
@@ -557,6 +569,7 @@ Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
     return kTRUE;
   }
   
+  //TODO Proper handle of tuneMConData for other dEdx sources
   
   Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used)
   Char_t ncl[3];        //same
@@ -572,8 +585,8 @@ Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
   if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows);
 
   //check if we cross a bad OROC in which case we reject
-  EChamberStatus trackStatus = TrackStatus(track,2);
-  if (trackStatus==kChamberOff || trackStatus==kChamberLowGain)
+  EChamberStatus trackOROCStatus = TrackStatus(track,2);
+  if (trackOROCStatus==kChamberOff || trackOROCStatus==kChamberLowGain)
   {
     return kFALSE;
   }
@@ -582,6 +595,7 @@ Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
   {
     case kdEdxOROC:
       {
+        if (trackOROCStatus==kChamberInvalid) return kFALSE; //never reached OROC
         dEdx = signal[3];
         nPoints = ncl[2]+ncl[1];
         gainScenario = kOROChigh;
@@ -852,63 +866,16 @@ Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0
 
 
 //_________________________________________________________________________
-Bool_t AliTPCPIDResponse::sectorNumbersInOut(const AliVTrack* track,
-                                             Double_t innerRadius,
-                                             Double_t outerRadius,
-                                             Float_t& inphi, 
+Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner,
+                                             Double_t* trackPositionOuter,
+                                             Float_t& inphi,
                                              Float_t& outphi,
-                                             Int_t& in, 
-                                             Int_t& out ) const
+                                             Int_t& in,
+                                                                                                                                                                                                                                     Int_t& out ) const
 {
   //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses
   //for OROC chamber numbers add 36
   //returned angles are between (0,2pi)
-  
-  Double_t trackPositionInner[3]; 
-  Double_t trackPositionOuter[3]; 
-
-  Bool_t trackAtInner=kTRUE;
-  Bool_t trackAtOuter=kTRUE;
-  const AliExternalTrackParam* ip = track->GetInnerParam();
-  
-  //if there is no inner param this could mean we're using the AOD track,
-  //we still can extrapolate from the vertex - so use those params.
-  if (ip) track=ip;
-
-  if (!track->GetXYZAt(innerRadius, fMagField, trackPositionInner))
-    trackAtInner=kFALSE;
-
-  if (!track->GetXYZAt(outerRadius, fMagField, trackPositionOuter))
-    trackAtOuter=kFALSE;
-
-  if (!trackAtInner)
-  {
-    //if we dont even enter inner radius we do nothing
-    inphi=0.0;
-    outphi=0.0;
-    in=0;
-    out=0;
-    return kFALSE;
-  }
-
-  if (!trackAtOuter)
-  {
-    //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
-    Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
-    Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
-    if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
-    {
-      //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
-    }
-    else
-    {
-      inphi=0.0;
-      outphi=0.0;
-      in=0;
-      out=0;
-      return kFALSE;
-    }
-  }
 
   inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]);
   outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]);
@@ -927,7 +894,8 @@ Bool_t AliTPCPIDResponse::sectorNumbersInOut(const AliVTrack* track,
   }
   return kTRUE;
 }
-
+    
+    
 //_____________________________________________________________________________
 Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
 {
@@ -955,14 +923,62 @@ AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack
   Float_t outphi=0.;
   Float_t innerRadius = (layer==1)?83.0:133.7;
   Float_t outerRadius = (layer==1)?133.5:247.7;
-  if (!sectorNumbersInOut(track, 
-                          innerRadius, 
-                          outerRadius, 
+
+  /////////////////////////////////////////////////////////////////////////////
+  //find out where track enters and leaves the layer.
+  //
+  Double_t trackPositionInner[3]; 
+  Double_t trackPositionOuter[3]; 
+  
+  //if there is no inner param this could mean we're using the AOD track,
+  //we still can extrapolate from the vertex - so use those params.
+  const AliExternalTrackParam* ip = track->GetInnerParam();
+  if (ip) track=ip;
+
+  Bool_t trackAtInner = track->GetXYZAt(innerRadius, fMagField, trackPositionInner);
+  Bool_t trackAtOuter = track->GetXYZAt(outerRadius, fMagField, trackPositionOuter);
+
+  if (!trackAtInner)
+  {
+    //if we dont even enter inner radius we do nothing and return invalid
+    inphi=0.0;
+    outphi=0.0;
+    in=0;
+    out=0;
+    return kChamberInvalid;
+  }
+
+  if (!trackAtOuter)
+  {
+    //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
+    Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
+    Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
+    if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
+    {
+      //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
+    }
+    else
+    {
+      inphi=0.0;
+      outphi=0.0;
+      in=0;
+      out=0;
+      return kChamberInvalid;
+    }
+  }
+
+
+  if (!sectorNumbersInOut(trackPositionInner, 
+                          trackPositionOuter, 
                           inphi, 
                           outphi, 
                           in, 
-                          out)) return kChamberOff;
+                          out)) return kChamberInvalid;
 
+  /////////////////////////////////////////////////////////////////////////////
+  //now we have the location of the track we can check 
+  //if it is in a good/bad chamber
+  //
   Bool_t sideA = kTRUE;
   
   if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE;
@@ -1000,6 +1016,7 @@ AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack
   Float_t trackLengthInBad = 0.;
   Float_t trackLengthInLowGain = 0.;
   Float_t trackLengthTotal = TMath::Abs(outphi-inphi);
+  Float_t lengthFractionInBadSectors = 0.;
 
   const Float_t sectorWidth = TMath::TwoPi()/18.;  
   
@@ -1016,12 +1033,21 @@ AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack
     else {deltaPhi=sectorWidth;}
     
     Float_t v = fVoltageMap[(layer==1)?(j):(j+36)];
-    if (v<=fBadIROCthreshhold) trackLengthInBad+=deltaPhi;
-    if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold) trackLengthInLowGain+=deltaPhi;
+    if (v<=fBadIROCthreshhold) 
+    { 
+      trackLengthInBad+=deltaPhi; 
+      lengthFractionInBadSectors=1.;
+    }
+    if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold) 
+    {
+      trackLengthInLowGain+=deltaPhi;
+      lengthFractionInBadSectors=1.;
+    }
   }
 
   //for now low gain and bad (off) chambers are treated equally
-  Float_t lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;
+  if (trackLengthTotal>0)
+    lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;
 
   //printf("### side: %s, pt: %.2f, pz: %.2f, in: %i, out: %i, phiIN: %.2f, phiOUT: %.2f, rIN: %.2f, rOUT: %.2f\n",(sideA)?"A":"C",track->Pt(),track->Pz(),in,out,inphi,outphi,innerRadius,outerRadius);
   
@@ -1091,4 +1117,4 @@ Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Do
   position[1]=b[1]+b[1]*TMath::Abs(r)/norm;
   position[2]=0.;
   return kTRUE;
-}
\ No newline at end of file
+}
index ce09635..a39ff6e 100644 (file)
@@ -126,8 +126,8 @@ public:
                                AliPID::EParticleType species,
                                ETPCdEdxSource dedxSource,
                                Double_t& dEdx, Int_t& nPoints, ETPCgainScenario& gainScenario, TSpline3** responseFunction) const;
-  Bool_t sectorNumbersInOut(const AliVTrack* track, 
-                            Double_t innerRadius, Double_t outerRadius, 
+  Bool_t sectorNumbersInOut(Double_t* trackPositionInner,
+                            Double_t* trackPositionOuter,
                             Float_t& phiIn, Float_t& phiOut, 
                             Int_t& in, Int_t& out ) const;
   AliTPCPIDResponse::EChamberStatus TrackStatus(const AliVTrack* track, Int_t layer) const;