3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright © ALICE HLT Group
6 #include "AliL3StandardIncludes.h"
8 #include "AliL3Logging.h"
9 #include "AliL3Modeller.h"
10 #include "AliL3MemHandler.h"
11 #include "AliL3TrackArray.h"
12 #include "AliL3ModelTrack.h"
13 #include "AliL3DigitData.h"
14 #include "AliL3Transform.h"
15 #include "AliL3SpacePointData.h"
18 #include "AliL3FileHandler.h"
25 //_____________________________________________________________
28 // Class for modeling TPC data.
30 // This performs the cluster finding, based on track parameters.
31 // Basically it propagates the tracks to all padrows, and looks
32 // for a corresponding cluster. For the moment only cog is calculated,
33 // and no deconvolution is done.
35 ClassImp(AliL3Modeller)
37 AliL3Modeller::AliL3Modeller()
46 SetMaxClusterRange(0,0);
51 AliL3Modeller::~AliL3Modeller()
61 void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,Bool_t houghtracks,Bool_t binary)
65 fHoughTracks=houghtracks;
67 sprintf(fPath,"%s",path);
69 fTracks = new AliL3TrackArray("AliL3ModelTrack");
72 AliL3MemHandler *file = new AliL3MemHandler();
74 sprintf(fname,"%s/tracks_tr_%d_0.raw",trackdata,fSlice); //output tracks from the tracker (no merging)
76 sprintf(fname,"%s/tracks_ho_%d.raw",trackdata,fSlice);
77 //sprintf(fname,"%s/tracks_ho_%d_%d.raw",trackdata,fSlice,fPatch);
78 if(!file->SetBinaryInput(fname))
80 cerr<<"AliL3Modeller::Init : Error opening trackfile: "<<fname<<endl;
83 file->Binary2TrackArray(fTracks);
84 file->CloseBinaryInput();
90 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
92 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
94 track->Init(fSlice,fPatch);
96 //Only if the tracks has been merged across sector boundaries:
98 //track->Rotate(fSlice,kTRUE); //!!!!!!!!!!!!!!!!!!!
100 track->CalculateHelix();
103 Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
104 Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
105 Int_t bounds = ntimes*npads;
106 fRow = new Digit[bounds];
110 AliL3DigitRowData *digits=0;
112 fMemHandler = new AliL3FileHandler();
113 fMemHandler->Init(slice,patch);
116 sprintf(fname,"%s/digitfile.root",fPath);
117 fMemHandler->SetAliInput(fname);
118 digits = fMemHandler->AliAltroDigits2Memory(ndigits);
122 sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
123 if(!fMemHandler->SetBinaryInput(fname))
125 cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
128 digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
131 fMemHandler = new AliL3MemHandler();
132 fMemHandler->Init(slice,patch);
135 cerr<<"AliL3Modeller::Init : Compile with AliROOT if you want rootfile as input"<<endl;
140 sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
141 if(!fMemHandler->SetBinaryInput(fname))
143 cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
147 digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
150 SetInputData(digits);
153 void AliL3Modeller::FindClusters()
156 cout<<"AliL3Modeller::FindClusters : Processing slice "<<fSlice<<" patch "<<fPatch<<endl;
159 cerr<<"AliL3Modeller::Process : No tracks"<<endl;
164 cerr<<"AliL3Modeller::Process : No data "<<endl;
168 AliL3DigitRowData *rowPt = fRowData;
169 AliL3DigitData *digPt=0;
174 ClusterRegion region[200];
176 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
178 if(i != (Int_t)rowPt->fRow)
180 cerr<<"AliL3Modeller::FindClusters : Mismatching rownumbering "<<i<<" "<<rowPt->fRow<<endl;
184 memset((void*)fRow,0,(AliL3Transform::GetNTimeBins()+1)*(AliL3Transform::GetNPads(i)+1)*sizeof(Digit));
185 digPt = (AliL3DigitData*)rowPt->fDigitData;
186 //cout<<"Loading row "<<i<<" with "<<(Int_t)rowPt->fNDigit<<" digits"<<endl;
187 for(UInt_t j=0; j<rowPt->fNDigit; j++)
190 time = digPt[j].fTime;
191 charge = digPt[j].fCharge;
192 fRow[(AliL3Transform::GetNTimeBins()+1)*pad + time].fCharge = charge;
193 fRow[(AliL3Transform::GetNTimeBins()+1)*pad + time].fUsed = kFALSE;
194 //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
197 for(Int_t k=0; k<fTracks->GetNTracks(); k++)
199 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
202 if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetNOverlaps(i)>0)//track->GetOverlap(i)>=0)
204 //cout<<"Track "<<k<<" is empty on row "<<i<<" "<<track->GetPadHit(i)<<" "<<track->GetTimeHit(i)<<endl;
205 track->SetCluster(i,0,0,0,0,0,0); //The track has left the patch, or it is overlapping
209 Int_t minpad,mintime,maxpad,maxtime;
210 minpad = mintime = 999;
211 maxpad = maxtime = 0;
213 memset(&cluster,0,sizeof(Cluster));
214 LocateCluster(track,region,minpad,maxpad);//,mintime,maxtime);
215 if(maxpad - minpad + 1 > fMaxPads || // maxtime - mintime + 1 > fMaxTimebins ||
216 maxpad - minpad < 1) // || maxtime - mintime < 1)
218 //cout<<"Cluster not found on row "<<i<<" maxpad "<<maxpad<<" minpad "<<minpad<<" maxtime "<<maxtime<<" mintime "<<mintime
219 // <<" padhit "<<track->GetPadHit(i)<<" timehit "<<track->GetTimeHit(i)<<endl;
221 track->SetCluster(i,0,0,0,0,0,0);
226 for(pad=minpad; pad<=maxpad; pad++)
229 for(time=region[pad].mintime; time<=region[pad].maxtime; time++)
231 charge = fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge;
232 if(!charge) continue;
233 if(fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed == kTRUE)
237 //Update the cluster parameters with this timebin
238 cluster.fTime += time*charge;
239 cluster.fPad += pad*charge;
240 cluster.fCharge += charge;
241 cluster.fSigmaY2 += pad*pad*charge;
242 cluster.fSigmaZ2 += time*time*charge;
243 fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed = kTRUE;
248 FillCluster(track,&cluster,i,npads);
251 fMemHandler->UpdateRowPointer(rowPt);
253 //cout<<"done processing"<<endl;
257 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
259 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
261 if(track->GetNClusters() != AliL3Transform::GetNRows(fPatch))
262 cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<AliL3Transform::GetNRows(fPatch)<<endl<<endl;
268 void AliL3Modeller::LocateCluster(AliL3ModelTrack *track,ClusterRegion *region,Int_t &padmin,Int_t &padmax)
270 //Set the cluster range
271 //This method searches for _all_ nonzeros timebins which are neigbours.
272 //This makes it rather impractical when dealing with high occupancy,
273 //because then you might have very large "cluster" areas from low
274 //pt electrons/noise.
276 Int_t row=fCurrentPadRow,charge,prtmin=0,prtmax=999;
277 Int_t hitpad = (Int_t)rint(track->GetPadHit(row));
278 Int_t hittime = (Int_t)rint(track->GetTimeHit(row));
279 Int_t tmin = hittime;
282 Int_t clustercharge=0;
285 Int_t npads=0,middlemax=tmax,middlemin=tmin;
289 Int_t time = hittime;
300 else if(pad >= AliL3Transform::GetNPads(row))
302 padmax = AliL3Transform::GetNPads(row)-1;
309 //cout<<"Starting to look in pad "<<pad<<" time "<<time<<endl;
317 else if(time >= AliL3Transform::GetNTimeBins())
319 //timemax = AliL3Transform::GetNTimeBins()-1;
322 charge = fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge;
324 //cout<<"charge "<<charge<<" at pad "<<pad<<" time "<<time<<endl;
327 clustercharge+=charge;
342 //if(abs(time - hittime) < fTimeSearch && padpr == kFALSE)//Keep looking
343 if(time > prtmin && npads!=0)
351 //else if(abs(time-hittime) < fTimeSearch && padpr == kFALSE)//Keep looking
352 else if(time < prtmax && npads != 0)
364 //cout<<"tmax "<<tmax<<" tmin "<<tmin<<" prtmin "<<prtmin<<" ptrmax "<<prtmax<<endl;
366 if(padpr && tmax >= prtmin && tmin <= prtmax)//Sequence is overlapping with the previous
369 //cout<<"Incrementing pad "<<endl;
372 region[pad].mintime=tmin;
373 region[pad].maxtime=tmax;
397 if(abs(pad-hitpad)<fPadSearch && clustercharge == 0)
402 //cout<<"Setting new pad "<<hitpad+1<<endl;
412 if(abs(pad-hitpad)<fPadSearch && clustercharge==0)
423 void AliL3Modeller::FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t row,Int_t npads)
425 if(cluster->fCharge==0)
427 track->SetCluster(row,0,0,0,0,0,0);
430 Float_t fcharge = (Float_t)cluster->fCharge;
431 Float_t fpad = ((Float_t)cluster->fPad/fcharge);
432 Float_t ftime = ((Float_t)cluster->fTime/fcharge);
433 Float_t sigmaY2,sigmaZ2;
434 CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
435 track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2,npads);
438 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
439 track->SetClusterLabel(row,trackID);
445 void AliL3Modeller::FillZeros(AliL3DigitRowData *rowPt,Bool_t reversesign)
447 //Fill zero where data has been used.
449 AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
450 for(UInt_t j=0; j<rowPt->fNDigit; j++)
452 Int_t pad = digPt[j].fPad;
453 Int_t time = digPt[j].fTime;
454 if(fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed==kTRUE)
458 if(digPt[j].fCharge < 1024)
459 digPt[j].fCharge += 1024;
462 digPt[j].fCharge = 0;
467 void AliL3Modeller::WriteRemaining()
469 //Write remaining (nonzero) digits to file.
471 AliL3DigitRowData *rowPt;
472 rowPt = (AliL3DigitRowData*)fRowData;
474 Int_t *ndigits=new Int_t[(AliL3Transform::GetNRows(fPatch))];
475 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
477 AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
478 ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]=0;
479 for(UInt_t j=0; j<rowPt->fNDigit; j++)
481 if(digPt[j].fCharge==0) continue;
483 ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]++;
485 //cout<<"Difference "<<(int)ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]<<" "<<(int)rowPt->fNDigit<<endl;
486 fMemHandler->UpdateRowPointer(rowPt);
489 Int_t size = digitcount*sizeof(AliL3DigitData) + AliL3Transform::GetNRows(fPatch)*sizeof(AliL3DigitRowData);
490 Byte_t *data = new Byte_t[size];
492 AliL3DigitRowData *tempPt = (AliL3DigitRowData*)data;
493 rowPt = (AliL3DigitRowData*)fRowData;
495 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
499 tempPt->fNDigit = ndigits[(i-AliL3Transform::GetFirstRow(fPatch))];
500 AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
501 for(UInt_t j=0; j<rowPt->fNDigit; j++)
503 if(digPt[j].fCharge==0) continue;
504 if(localcount >= ndigits[(i-AliL3Transform::GetFirstRow(fPatch))])
506 cerr<<"AliL3Modeller::WriteRemaining : Digitarray out of range!!"<<endl;
509 tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
510 tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
511 tempPt->fDigitData[localcount].fTime = digPt[j].fTime;
515 if(ndigits[(i-AliL3Transform::GetFirstRow(fPatch))] != localcount)
517 cerr<<"AliL3Modeller::WriteRemaining : Mismatch in digitcount"<<endl;
520 fMemHandler->UpdateRowPointer(rowPt);
521 Byte_t *tmp = (Byte_t*)tempPt;
522 Int_t size = sizeof(AliL3DigitRowData) + ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]*sizeof(AliL3DigitData);
524 tempPt = (AliL3DigitRowData*)tmp;
528 AliL3MemHandler *mem = new AliL3MemHandler();
529 sprintf(fname,"%s/comp/remains_%d_%d.raw",fPath,fSlice,fPatch);
530 mem->Init(fSlice,fPatch);
531 mem->SetBinaryOutput(fname);
532 mem->Memory2CompBinary((UInt_t)AliL3Transform::GetNRows(fPatch),(AliL3DigitRowData*)data);
533 mem->CloseBinaryOutput();
539 void AliL3Modeller::RemoveBadTracks()
541 //Remove tracsk which should not be included in the compression scheme.
543 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
545 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
548 if(track->GetPt() < 0.08)
555 if(track->GetNHits() < fTrackThreshold)
562 void AliL3Modeller::CalculateCrossingPoints()
565 cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
568 cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
574 for(Int_t i=AliL3Transform::GetLastRow(fPatch); i>=AliL3Transform::GetFirstRow(fPatch); i--)
576 for(Int_t j=0; j<fTracks->GetNTracks(); j++)
578 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
581 if(!track->GetCrossingPoint(i,hit))
583 //cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
584 // " pt "<<track->GetPt()<<
585 // " tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
586 //fTracks->Remove(j);
587 track->SetPadHit(i,-1);
588 track->SetTimeHit(i,-1);
591 //cout<<"X "<<hit[0]<<" Y "<<hit[1]<<" Z "<<hit[2]<<" tgl "<<track->GetTgl()<<endl;
593 AliL3Transform::Slice2Sector(fSlice,i,sector,row);
594 AliL3Transform::Local2Raw(hit,sector,row);
595 //cout<<"Pad "<<hit[1]<<" time "<<hit[2]<<" in sector "<<sector<<" row "<<row<<endl;
596 if(hit[1]<0 || hit[1]>AliL3Transform::GetNPads(i) ||
597 hit[2]<0 || hit[2]>AliL3Transform::GetNTimeBins())
598 {//Track is leaving the patch, so flag the track hits (<0)
599 track->SetPadHit(i,-1);
600 track->SetTimeHit(i,-1);
604 track->SetPadHit(i,hit[1]);
605 track->SetTimeHit(i,hit[2]);
606 track->CalculateClusterWidths(i);
608 Double_t beta = track->GetCrossingAngle(i);
609 track->SetCrossingAngleLUT(i,beta);
611 //if(hit[1]<0 || hit[2]>445)
612 //if(hit[2]<0 || hit[2]>445)
613 //cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
614 //cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
619 cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
622 void AliL3Modeller::CheckForOverlaps(Float_t dangle,Int_t *rowrange)
624 //Flag the tracks that overlap
627 cout<<"Checking for overlaps on "<<fTracks->GetNTracks()<<endl;
630 for(Int_t k=AliL3Transform::GetFirstRow(fPatch); k<=AliL3Transform::GetLastRow(fPatch); k++)
634 if(k < rowrange[0]) continue;
635 if(k > rowrange[1]) break;
637 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
639 AliL3ModelTrack *track1 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
640 if(!track1) continue;
641 if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0) continue;
643 for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
645 AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
646 if(!track2) continue;
647 if(track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0) continue;
649 if(abs((Int_t)rint(track1->GetPadHit(k))-(Int_t)rint(track2->GetPadHit(k))) <= fPadOverlap &&
650 abs((Int_t)rint(track1->GetTimeHit(k))-(Int_t)rint(track2->GetTimeHit(k))) <= fTimeOverlap)
652 if(dangle>0 && fabs(track1->GetCrossingAngleLUT(k) - track2->GetCrossingAngleLUT(k)) < dangle)
655 //cout<<"row "<<k<<" "<<i<<" "<<j<<" "<<track1->GetPadHit(k)<<" "<<track2->GetPadHit(k)<<" "<<fabs(track1->GetCrossingAngleLUT(k) - track2->GetCrossingAngleLUT(k))<<endl;
658 track1->SetOverlap(k,j);
666 cout<<"and there are "<<fTracks->GetNTracks()<<" track left"<<endl;
667 //cout<<"found "<<counter<<" done"<<endl;
671 void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
676 padw = AliL3Transform::GetPadPitchWidth(fPatch);
678 Float_t charge = (Float_t)cl->fCharge;
679 Float_t pad = (Float_t)cl->fPad/charge;
680 Float_t time = (Float_t)cl->fTime/charge;
681 Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
683 //Save the sigmas in pad and time:
685 sigmaY2 = (s2);// + 1./12);//*padw*padw;
687 /*Constants added by offline
690 sigmaY2 = sigmaY2*0.108;
692 sigmaY2 = sigmaY2*2.07;
696 s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
697 timew = AliL3Transform::GetZWidth();
698 sigmaZ2 = (s2);// +1./12);//*timew*timew;
703 Constants added by offline
706 sigmaZ2 = sigmaZ2*0.169;
708 sigmaZ2 = sigmaZ2*1.77;
713 void AliL3Modeller::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
716 AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fRowData;
718 trackID[0]=trackID[1]=trackID[2]=-2;
720 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
722 if(rowPt->fRow < (UInt_t)fCurrentPadRow)
724 AliL3MemHandler::UpdateRowPointer(rowPt);
727 AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
728 for(UInt_t j=0; j<rowPt->fNDigit; j++)
730 Int_t cpad = digPt[j].fPad;
731 Int_t ctime = digPt[j].fTime;
732 if(cpad != pad) continue;
733 if(ctime != time) continue;
734 //if(cpad != pad && ctime != ctime) continue;
735 //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
736 trackID[0] = digPt[j].fTrackID[0];
737 trackID[1] = digPt[j].fTrackID[1];
738 trackID[2] = digPt[j].fTrackID[2];
740 //cout<<"Reading trackID "<<trackID[0]<<endl;