+//______________________________________________________
+
+Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const
+{
+ //
+ // Returns number of the innermost SENSITIVE propagation layer
+ //
+
+ return GetLayerNumber(0);
+}
+
+//______________________________________________________
+
+Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const
+{
+ //
+ // Returns number of the outermost SENSITIVE time bin
+ //
+
+ return GetLayerNumber(GetNumberOfTimeBins() - 1);
+}
+
+//______________________________________________________
+
+Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const
+{
+ //
+ // Returns number of SENSITIVE time bins
+ //
+
+ Int_t tb, layer;
+ for(tb = kMaxTimeBinIndex-1; tb >=0; tb--) {
+ layer = GetLayerNumber(tb);
+ if(layer>=0) break;
+ }
+ return tb+1;
+}
+
+//______________________________________________________
+
+void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
+{
+ //
+ // Insert layer <pl> in fLayers array.
+ // Layers are sorted according to X coordinate.
+
+ if ( fN == ((Int_t) kMaxLayersPerSector)) {
+ printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
+ return;
+ }
+ if (fN==0) {fLayers[fN++] = pl; return;}
+ Int_t i=Find(pl->GetX());
+
+ memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
+ fLayers[i]=pl; fN++;
+
+}
+
+//______________________________________________________
+
+Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const
+{
+ //
+ // Returns index of the propagation layer nearest to X
+ //
+
+ if (x <= fLayers[0]->GetX()) return 0;
+ if (x > fLayers[fN-1]->GetX()) return fN;
+ Int_t b=0, e=fN-1, m=(b+e)/2;
+ for (; b<e; m=(b+e)/2) {
+ if (x > fLayers[m]->GetX()) b=m+1;
+ else e=m;
+ }
+ return m;
+}
+
+//______________________________________________________
+void AliTRDtracker::AliTRDpropagationLayer::SetZ(Double_t* center, Double_t *w, Double_t *wsensitive )
+{
+ //
+ // set centers and the width of sectors
+ for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
+ fZc[icham] = center[icham];
+ fZmax[icham] = w[icham];
+ fZmaxSensitive[icham] = wsensitive[icham];
+ // printf("chamber\t%d\tzc\t%f\tzmax\t%f\tzsens\t%f\n",icham,fZc[icham],fZmax[icham],fZmaxSensitive[icham]);
+ }
+}
+//______________________________________________________
+
+void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
+{
+ //
+ // set centers and the width of sectors
+ fHole = kFALSE;
+ for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
+ fIsHole[icham] = holes[icham];
+ if (holes[icham]) fHole = kTRUE;
+ }
+}
+
+
+
+void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
+ Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &radLength,
+ Bool_t &lookForCluster) const
+{
+ //
+ // Returns radial step <dx>, density <rho>, rad. length <radLength>,
+ // and sensitivity <lookForCluster> in point <y,z>
+ //
+
+ Double_t alpha = AliTRDgeometry::GetAlpha();
+ Double_t ymax = fX*TMath::Tan(0.5*alpha);
+
+
+ dx = fdX;
+ rho = fRho;
+ radLength = fX0;
+ lookForCluster = kFALSE;
+ Bool_t cross =kFALSE;
+ //
+ //
+ if ( (ymax-TMath::Abs(y))<3.){ //cross material
+ rho*=40.;
+ radLength*=40.;
+ cross=kTRUE;
+ }
+ //
+ // check dead regions in sensitive volume
+ //
+ Int_t zone=-1;
+ for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
+ if (TMath::Abs(z - fZc[ch]) > fZmax[ch]) continue; //not in given zone
+ //
+ if (TMath::Abs(z - fZc[ch]) < fZmaxSensitive[ch]){
+ if (fTimeBinIndex>=0) lookForCluster = !(fIsHole[zone]);
+ if(TMath::Abs(y) > fYmaxSensitive){
+ lookForCluster = kFALSE;
+ }
+ if (fIsHole[zone]) {
+ //if hole
+ rho = 1.29e-3;
+ radLength = 36.66;
+ }
+ }else{
+ rho = 2.7; radLength = 24.01; //aluminium in between
+ }
+ }
+ //
+ if (fTimeBinIndex>=0) return;
+ //
+ //
+ // check hole
+ if (fHole==kFALSE) return;
+ //
+ for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
+ if (TMath::Abs(z - fZc[ch]) < fZmax[ch]){
+ if (fIsHole[ch]) {
+ //if hole
+ rho = 1.29e-3;
+ radLength = 36.66;
+ }
+ }
+ }
+ return;
+}
+
+Int_t AliTRDtracker::AliTRDpropagationLayer::GetZone( Double_t z) const
+{
+ //
+ //
+ if (fTimeBinIndex < 0) return -20; //unknown
+ Int_t zone=-10; // dead zone
+ for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
+ if(TMath::Abs(z - fZc[ch]) < fZmax[ch])
+ zone = ch;
+ }
+ return zone;
+}
+
+
+//______________________________________________________
+
+void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c,
+ UInt_t index) {
+
+// Insert cluster in cluster array.
+// Clusters are sorted according to Y coordinate.
+
+ if(fTimeBinIndex < 0) {
+ printf("*** attempt to insert cluster into non-sensitive time bin!\n");
+ return;
+ }
+
+ if (fN== (Int_t) kMaxClusterPerTimeBin) {
+ printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n");
+ return;
+ }
+ if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
+ Int_t i=Find(c->GetY());
+ memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
+ memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
+ fIndex[i]=index; fClusters[i]=c; fN++;
+}
+
+//______________________________________________________
+
+Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Double_t y) const {
+
+// Returns index of the cluster nearest in Y
+
+ if (y <= fClusters[0]->GetY()) return 0;
+ if (y > fClusters[fN-1]->GetY()) return fN;
+ Int_t b=0, e=fN-1, m=(b+e)/2;
+ for (; b<e; m=(b+e)/2) {
+ if (y > fClusters[m]->GetY()) b=m+1;
+ else e=m;
+ }
+ return m;
+}
+
+//---------------------------------------------------------
+
+Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
+//
+// Returns correction factor for tilted pads geometry
+//
+
+ AliTRDpadPlane *padPlane = fPar->GetPadPlane(0,0);
+ Double_t h01 = sin(TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
+ Int_t det = c->GetDetector();
+ Int_t plane = fGeom->GetPlane(det);
+
+ //if((plane == 1) || (plane == 3) || (plane == 5)) h01=-h01;
+ if((plane == 0) || (plane == 2) || (plane == 4)) h01=-h01;
+
+ if(fNoTilt) h01 = 0;
+
+ return h01;
+}
+
+
+void AliTRDtracker::CookdEdxTimBin(AliTRDtrack& TRDtrack)
+{
+ // *** ADDED TO GET MORE INFORMATION FOR TRD PID ---- PS
+ // 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)
+
+ // const Int_t kNPlane = AliTRDgeometry::Nplan();
+ // const Int_t kNPlane = 6;
+ Double_t clscharge[kNPlane], maxclscharge[kNPlane];
+ Int_t nCluster[kNPlane], timebin[kNPlane];
+
+ //Initialization of cluster charge per plane.
+ for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
+ clscharge[iPlane] = 0.0;
+ nCluster[iPlane] = 0;
+ timebin[iPlane] = -1;
+ maxclscharge[iPlane] = 0.0;
+ }
+
+ // Loop through all clusters associated to track TRDtrack
+ 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);
+ AliTRDcluster *TRDcluster = (AliTRDcluster *) GetCluster(index);
+ if (!TRDcluster) continue;
+ Int_t tb = TRDcluster->GetLocalTimeBin();
+ if (!tb) continue;
+ Int_t detector = TRDcluster->GetDetector();
+ Int_t iPlane = fGeom->GetPlane(detector);
+ clscharge[iPlane] = clscharge[iPlane]+charge;
+ if(charge > maxclscharge[iPlane]) {
+ maxclscharge[iPlane] = charge;
+ timebin[iPlane] = tb;
+ }
+ nCluster[iPlane]++;
+ } // end of loop over cluster
+
+ // Setting the fdEdxPlane and fTimBinPlane variabales
+ Double_t Total_ch = 0;
+ for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
+ // Quality control of TRD track.
+ if (nCluster[iPlane]<= 5) {
+ clscharge[iPlane]=0.0;
+ timebin[iPlane]=-1;
+ }
+ if (nCluster[iPlane]) clscharge[iPlane] /= nCluster[iPlane];
+ TRDtrack.SetPIDsignals(clscharge[iPlane], iPlane);
+ TRDtrack.SetPIDTimBin(timebin[iPlane], iPlane);
+ Total_ch= Total_ch+clscharge[iPlane];
+ }
+ // Int_t i;
+ // Int_t nc=TRDtrack.GetNumberOfClusters();
+ // Float_t dedx=0;
+ // for (i=0; i<nc; i++) dedx += TRDtrack.GetClusterdQdl(i);
+ // dedx /= nc;
+ // for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
+ // TRDtrack.SetPIDsignals(dedx, iPlane);
+ // TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
+ // }
+
+} // end of function
+
+
+Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters)
+{
+ //
+ //
+ // try to find nearest clusters to the track in timebins from t0 to t1
+ //
+ //
+ Double_t x[1000],yt[1000],zt[10000];
+ Double_t dz[1000],dy[10000];
+ Int_t indexes[2][10000];
+ AliTRDcluster *cl[2][10000];
+
+ for (Int_t it=t0;it<t1; it++){
+ clusters[it]=-2;
+ indexes[0][it]=-2; //reset indexes1
+ indexes[1][it]=-2; //reset indexes1
+ x[it]=0;
+ yt[it]=0;
+ zt[it]=0;
+ dz[it]=0;
+ dy[it]=0;
+ cl[0][it]=0;
+ cl[1][it]=0;
+ }
+ //
+ Double_t x0 = track->GetX();
+ Double_t sy2=ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
+ Double_t sz2=ExpectedSigmaZ2(x0,track->GetTgl());
+ Double_t road = 10.*sqrt(track->GetSigmaY2() + sy2);
+ Int_t nall=0;
+ Int_t nfound=0;
+
+ for (Int_t it=t0;it<t1;it++){
+ Double_t maxChi2=fgkMaxChi2;
+ AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(it));
+ if (timeBin==0) continue; // no indexes1
+ Int_t maxn = timeBin;
+ x[it] = timeBin.GetX();
+ Double_t y,z;
+ if (!track->GetProlongation(x[it],y,z)) continue;
+ yt[it]=y;
+ zt[it]=z;
+ Double_t chi2 =1000000;
+ nall++;
+ //
+ // find nearest cluster at given pad
+ for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
+ AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
+ Double_t h01 = GetTiltFactor(c);
+ if (c->GetY() > y+road) break;
+ if (c->IsUsed() > 0) continue;
+ if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
+ chi2=track->GetPredictedChi2(c,h01);
+ if (chi2 > maxChi2) continue;
+ maxChi2=chi2;
+ cl[0][it]=c;
+ indexes[0][it] =timeBin.GetIndex(i);
+ }
+ //
+ // find nearest cluster also in adjacent 2 pads
+ //
+ for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
+ AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
+ Double_t h01 = GetTiltFactor(c);
+ if (c->GetY() > y+road) break;
+ if (c->IsUsed() > 0) continue;
+ if((c->GetZ()-z)*(c->GetZ()-z) > 12 * sz2) continue;
+ chi2=track->GetPredictedChi2(c,h01);
+ if (chi2 > maxChi2) continue;
+ maxChi2=chi2;
+ cl[1][it]=c;
+ indexes[1][it]=timeBin.GetIndex(i);
+ if (!cl[0][it]) {
+ cl[0][it]=c;
+ indexes[0][it]=timeBin.GetIndex(i);
+ }
+ }
+ if (cl[0][it]||cl[1][it]) nfound++;
+ }
+ //
+ if (nfound<5) return 0;
+ //
+ // choose one of the variants
+ //
+ Int_t changes[2]={0,0};
+ Float_t sigmay[2]={1000,1000};
+ Float_t meany[2] ={1000,1000};
+ Float_t meanz[2] ={1000,1000};
+ Int_t sumall[2] ={0,0};
+ Int_t ngood[2] ={0,0};
+ Int_t nbad[2] ={0,0};
+ //
+ //
+ for (Int_t ih=0; ih<2;ih++){
+ Float_t lastz =-10000;
+ Float_t sumz =0;
+ Float_t sum =0;
+ Double_t sumdy = 0;
+ Double_t sumdy2= 0;
+ //
+ // how many changes ++ mean z ++mean y ++ sigma y
+ for (Int_t it=t0;it<t1;it++){
+ if (!cl[ih][it]) continue;
+ sumall[ih]++;
+ if (lastz<-9999) lastz = cl[ih][it]->GetZ();
+ if (TMath::Abs(lastz-cl[ih][it]->GetZ())>1) {
+ lastz = cl[ih][it]->GetZ();
+ changes[ih]++;
+ }
+ sumz = cl[ih][it]->GetZ();
+ sum++;
+ Double_t h01 = GetTiltFactor(cl[ih][it]);
+ dz[it] = cl[ih][it]->GetZ()- zt[it];
+ dy[it] = cl[ih][it]->GetY()+ dz[it]*h01 -yt[it];
+ sum++;
+ sumdy += dy[it];
+ sumdy2+= dy[it]*dy[it];
+ Int_t label = TMath::Abs(track->GetLabel());
+ if (cl[ih][it]->GetLabel(0)==label || cl[ih][it]->GetLabel(1)==label||cl[ih][it]->GetLabel(2)==label){
+ ngood[ih]++;
+ }
+ else{
+ nbad[ih]++;
+ }
+ }
+ if (sumall[ih]>4){
+ meanz[ih] = sumz/sum;
+ meany[ih] = sumdy/sum;
+ sigmay[ih] = TMath::Sqrt(sumdy2/sum-meany[ih]*meany[ih]);
+ }
+ }
+ Int_t best =0;
+ /*
+ if (sumall[0]<sumall[1]-2&&sigmay[1]<0.1){
+ if (sigmay[1]<sigmay[0]) best = 1; // if sigma is better
+ }
+ */
+ Float_t quality0 = (sigmay[0]+TMath::Abs(meany[0]))*(1.+Float_t(changes[0])/Float_t(sumall[0]));
+ Float_t quality1 = (sigmay[1]+TMath::Abs(meany[1]))*(1.+Float_t(changes[1])/Float_t(sumall[1]));
+
+ if (quality0>quality1){
+ best = 1;
+ }
+
+
+ //
+ for (Int_t it=t0;it<t1;it++){
+ if (!cl[best][it]) continue;
+ Double_t h01 = GetTiltFactor(cl[best][it]);
+ dz[it] = cl[best][it]->GetZ()- zt[it];
+ dy[it] = cl[best][it]->GetY()+ dz[it]*h01 -yt[it];
+ //
+ if (TMath::Abs(dy[it])<2.5*sigmay[best])
+ clusters[it] = indexes[best][it];
+ }
+
+ if (nbad[0]>4){
+ nbad[0] = nfound;
+ }
+ return nfound;
+}
+
+