3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright © ALICE HLT Group
5 //_____________________________________________________________
8 // Class for modeling TPC data.
10 // This performs the cluster finding, based on track parameters.
11 // Basically it propagates the tracks to all padrows, and looks
12 // for a corresponding cluster. For the moment only cog is calculated,
13 // and no deconvolution is done.
16 #include "AliL3StandardIncludes.h"
18 #include "AliL3Logging.h"
19 #include "AliL3Modeller.h"
20 #include "AliL3MemHandler.h"
21 #include "AliL3TrackArray.h"
22 #include "AliL3ModelTrack.h"
23 #include "AliL3DigitData.h"
24 #include "AliL3Transform.h"
25 #include "AliL3SpacePointData.h"
28 #include "AliL3FileHandler.h"
35 ClassImp(AliL3Modeller)
37 AliL3Modeller::AliL3Modeller()
39 // default constructor
47 SetMaxClusterRange(0,0);
52 AliL3Modeller::~AliL3Modeller()
63 void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,Bool_t houghtracks,Bool_t binary)
68 fHoughTracks=houghtracks;
70 sprintf(fPath,"%s",path);
72 fTracks = new AliL3TrackArray("AliL3ModelTrack");
75 AliL3MemHandler *file = new AliL3MemHandler();
77 sprintf(fname,"%s/tracks_tr_%d_0.raw",trackdata,fSlice); //output tracks from the tracker (no merging)
79 sprintf(fname,"%s/tracks_ho_%d.raw",trackdata,fSlice);
80 //sprintf(fname,"%s/tracks_ho_%d_%d.raw",trackdata,fSlice,fPatch);
81 if(!file->SetBinaryInput(fname))
83 cerr<<"AliL3Modeller::Init : Error opening trackfile: "<<fname<<endl;
86 file->Binary2TrackArray(fTracks);
87 file->CloseBinaryInput();
93 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
95 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
97 track->Init(fSlice,fPatch);
99 //Only if the tracks has been merged across sector boundaries:
101 //track->Rotate(fSlice,kTRUE); //!!!!!!!!!!!!!!!!!!!
103 track->CalculateHelix();
106 Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
107 Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
108 Int_t bounds = ntimes*npads;
109 fRow = new Digit[bounds];
113 AliL3DigitRowData *digits=0;
115 fMemHandler = new AliL3FileHandler();
116 fMemHandler->Init(slice,patch);
119 sprintf(fname,"%s/digitfile.root",fPath);
120 fMemHandler->SetAliInput(fname);
121 digits = fMemHandler->AliAltroDigits2Memory(ndigits);
125 sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
126 if(!fMemHandler->SetBinaryInput(fname))
128 cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
131 digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
134 fMemHandler = new AliL3MemHandler();
135 fMemHandler->Init(slice,patch);
138 cerr<<"AliL3Modeller::Init : Compile with AliROOT if you want rootfile as input"<<endl;
143 sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
144 if(!fMemHandler->SetBinaryInput(fname))
146 cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
150 digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
153 SetInputData(digits);
156 void AliL3Modeller::FindClusters()
160 cout<<"AliL3Modeller::FindClusters : Processing slice "<<fSlice<<" patch "<<fPatch<<endl;
163 cerr<<"AliL3Modeller::Process : No tracks"<<endl;
168 cerr<<"AliL3Modeller::Process : No data "<<endl;
172 AliL3DigitRowData *rowPt = fRowData;
173 AliL3DigitData *digPt=0;
178 ClusterRegion region[200];
180 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
182 if(i != (Int_t)rowPt->fRow)
184 cerr<<"AliL3Modeller::FindClusters : Mismatching rownumbering "<<i<<" "<<rowPt->fRow<<endl;
188 memset((void*)fRow,0,(AliL3Transform::GetNTimeBins()+1)*(AliL3Transform::GetNPads(i)+1)*sizeof(Digit));
189 digPt = (AliL3DigitData*)rowPt->fDigitData;
190 //cout<<"Loading row "<<i<<" with "<<(Int_t)rowPt->fNDigit<<" digits"<<endl;
191 for(UInt_t j=0; j<rowPt->fNDigit; j++)
194 time = digPt[j].fTime;
195 charge = digPt[j].fCharge;
196 fRow[(AliL3Transform::GetNTimeBins()+1)*pad + time].fCharge = charge;
197 fRow[(AliL3Transform::GetNTimeBins()+1)*pad + time].fUsed = kFALSE;
198 //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
201 for(Int_t k=0; k<fTracks->GetNTracks(); k++)
203 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
206 if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetNOverlaps(i)>0)//track->GetOverlap(i)>=0)
208 //cout<<"Track "<<k<<" is empty on row "<<i<<" "<<track->GetPadHit(i)<<" "<<track->GetTimeHit(i)<<endl;
209 track->SetCluster(i,0,0,0,0,0,0); //The track has left the patch, or it is overlapping
213 Int_t minpad,mintime,maxpad,maxtime;
214 minpad = mintime = 999;
215 maxpad = maxtime = 0;
217 memset(&cluster,0,sizeof(Cluster));
218 LocateCluster(track,region,minpad,maxpad);//,mintime,maxtime);
219 if(maxpad - minpad + 1 > fMaxPads || // maxtime - mintime + 1 > fMaxTimebins ||
220 maxpad - minpad < 1) // || maxtime - mintime < 1)
222 //cout<<"Cluster not found on row "<<i<<" maxpad "<<maxpad<<" minpad "<<minpad<<" maxtime "<<maxtime<<" mintime "<<mintime
223 // <<" padhit "<<track->GetPadHit(i)<<" timehit "<<track->GetTimeHit(i)<<endl;
225 track->SetCluster(i,0,0,0,0,0,0);
230 for(pad=minpad; pad<=maxpad; pad++)
233 for(time=region[pad].fMintime; time<=region[pad].fMaxtime; time++)
235 charge = fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge;
236 if(!charge) continue;
237 if(fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed == kTRUE)
241 //Update the cluster parameters with this timebin
242 cluster.fTime += time*charge;
243 cluster.fPad += pad*charge;
244 cluster.fCharge += charge;
245 cluster.fSigmaY2 += pad*pad*charge;
246 cluster.fSigmaZ2 += time*time*charge;
247 fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed = kTRUE;
252 FillCluster(track,&cluster,i,npads);
255 fMemHandler->UpdateRowPointer(rowPt);
257 //cout<<"done processing"<<endl;
261 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
263 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
265 if(track->GetNClusters() != AliL3Transform::GetNRows(fPatch))
266 cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<AliL3Transform::GetNRows(fPatch)<<endl<<endl;
272 void AliL3Modeller::LocateCluster(AliL3ModelTrack *track,ClusterRegion *region,Int_t &padmin,Int_t &padmax)
274 //Set the cluster range
275 //This method searches for _all_ nonzeros timebins which are neigbours.
276 //This makes it rather impractical when dealing with high occupancy,
277 //because then you might have very large "cluster" areas from low
278 //pt electrons/noise.
280 Int_t row=fCurrentPadRow,charge,prtmin=0,prtmax=999;
281 Int_t hitpad = (Int_t)rint(track->GetPadHit(row));
282 Int_t hittime = (Int_t)rint(track->GetTimeHit(row));
283 Int_t tmin = hittime;
286 Int_t clustercharge=0;
289 Int_t npads=0,middlemax=tmax,middlemin=tmin;
293 Int_t time = hittime;
304 else if(pad >= AliL3Transform::GetNPads(row))
306 padmax = AliL3Transform::GetNPads(row)-1;
313 //cout<<"Starting to look in pad "<<pad<<" time "<<time<<endl;
321 else if(time >= AliL3Transform::GetNTimeBins())
323 //timemax = AliL3Transform::GetNTimeBins()-1;
326 charge = fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge;
328 //cout<<"charge "<<charge<<" at pad "<<pad<<" time "<<time<<endl;
331 clustercharge+=charge;
346 //if(abs(time - hittime) < fTimeSearch && padpr == kFALSE)//Keep looking
347 if(time > prtmin && npads!=0)
355 //else if(abs(time-hittime) < fTimeSearch && padpr == kFALSE)//Keep looking
356 else if(time < prtmax && npads != 0)
368 //cout<<"tmax "<<tmax<<" tmin "<<tmin<<" prtmin "<<prtmin<<" ptrmax "<<prtmax<<endl;
370 if(padpr && tmax >= prtmin && tmin <= prtmax)//Sequence is overlapping with the previous
373 //cout<<"Incrementing pad "<<endl;
376 region[pad].fMintime=tmin;
377 region[pad].fMaxtime=tmax;
401 if(abs(pad-hitpad)<fPadSearch && clustercharge == 0)
406 //cout<<"Setting new pad "<<hitpad+1<<endl;
416 if(abs(pad-hitpad)<fPadSearch && clustercharge==0)
427 void AliL3Modeller::FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t row,Int_t npads)
430 if(cluster->fCharge==0)
432 track->SetCluster(row,0,0,0,0,0,0);
435 Float_t fcharge = (Float_t)cluster->fCharge;
436 Float_t fpad = ((Float_t)cluster->fPad/fcharge);
437 Float_t ftime = ((Float_t)cluster->fTime/fcharge);
438 Float_t sigmaY2,sigmaZ2;
439 CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
440 track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2,npads);
443 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
444 track->SetClusterLabel(row,trackID);
450 void AliL3Modeller::FillZeros(AliL3DigitRowData *rowPt,Bool_t reversesign)
452 //Fill zero where data has been used.
454 AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
455 for(UInt_t j=0; j<rowPt->fNDigit; j++)
457 Int_t pad = digPt[j].fPad;
458 Int_t time = digPt[j].fTime;
459 if(fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed==kTRUE)
463 if(digPt[j].fCharge < 1024)
464 digPt[j].fCharge += 1024;
467 digPt[j].fCharge = 0;
472 void AliL3Modeller::WriteRemaining()
474 //Write remaining (nonzero) digits to file.
476 AliL3DigitRowData *rowPt;
477 rowPt = (AliL3DigitRowData*)fRowData;
479 Int_t *ndigits=new Int_t[(AliL3Transform::GetNRows(fPatch))];
480 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
482 AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
483 ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]=0;
484 for(UInt_t j=0; j<rowPt->fNDigit; j++)
486 if(digPt[j].fCharge==0) continue;
488 ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]++;
490 //cout<<"Difference "<<(int)ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]<<" "<<(int)rowPt->fNDigit<<endl;
491 fMemHandler->UpdateRowPointer(rowPt);
494 Int_t size = digitcount*sizeof(AliL3DigitData) + AliL3Transform::GetNRows(fPatch)*sizeof(AliL3DigitRowData);
495 Byte_t *data = new Byte_t[size];
497 AliL3DigitRowData *tempPt = (AliL3DigitRowData*)data;
498 rowPt = (AliL3DigitRowData*)fRowData;
500 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
504 tempPt->fNDigit = ndigits[(i-AliL3Transform::GetFirstRow(fPatch))];
505 AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
506 for(UInt_t j=0; j<rowPt->fNDigit; j++)
508 if(digPt[j].fCharge==0) continue;
509 if(localcount >= ndigits[(i-AliL3Transform::GetFirstRow(fPatch))])
511 cerr<<"AliL3Modeller::WriteRemaining : Digitarray out of range!!"<<endl;
514 tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
515 tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
516 tempPt->fDigitData[localcount].fTime = digPt[j].fTime;
520 if(ndigits[(i-AliL3Transform::GetFirstRow(fPatch))] != localcount)
522 cerr<<"AliL3Modeller::WriteRemaining : Mismatch in digitcount"<<endl;
525 fMemHandler->UpdateRowPointer(rowPt);
526 Byte_t *tmp = (Byte_t*)tempPt;
527 Int_t size = sizeof(AliL3DigitRowData) + ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]*sizeof(AliL3DigitData);
529 tempPt = (AliL3DigitRowData*)tmp;
533 AliL3MemHandler *mem = new AliL3MemHandler();
534 sprintf(fname,"%s/comp/remains_%d_%d.raw",fPath,fSlice,fPatch);
535 mem->Init(fSlice,fPatch);
536 mem->SetBinaryOutput(fname);
537 mem->Memory2CompBinary((UInt_t)AliL3Transform::GetNRows(fPatch),(AliL3DigitRowData*)data);
538 mem->CloseBinaryOutput();
544 void AliL3Modeller::RemoveBadTracks()
546 //Remove tracsk which should not be included in the compression scheme.
548 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
550 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
553 if(track->GetPt() < 0.08)
560 if(track->GetNHits() < fTrackThreshold)
567 void AliL3Modeller::CalculateCrossingPoints()
569 // calculates crossing points
571 cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
574 cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
580 for(Int_t i=AliL3Transform::GetLastRow(fPatch); i>=AliL3Transform::GetFirstRow(fPatch); i--)
582 for(Int_t j=0; j<fTracks->GetNTracks(); j++)
584 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
587 if(!track->GetCrossingPoint(i,hit))
589 //cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
590 // " pt "<<track->GetPt()<<
591 // " tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
592 //fTracks->Remove(j);
593 track->SetPadHit(i,-1);
594 track->SetTimeHit(i,-1);
597 //cout<<"X "<<hit[0]<<" Y "<<hit[1]<<" Z "<<hit[2]<<" tgl "<<track->GetTgl()<<endl;
599 AliL3Transform::Slice2Sector(fSlice,i,sector,row);
600 AliL3Transform::Local2Raw(hit,sector,row);
601 //cout<<"Pad "<<hit[1]<<" time "<<hit[2]<<" in sector "<<sector<<" row "<<row<<endl;
602 if(hit[1]<0 || hit[1]>AliL3Transform::GetNPads(i) ||
603 hit[2]<0 || hit[2]>AliL3Transform::GetNTimeBins())
604 {//Track is leaving the patch, so flag the track hits (<0)
605 track->SetPadHit(i,-1);
606 track->SetTimeHit(i,-1);
610 track->SetPadHit(i,hit[1]);
611 track->SetTimeHit(i,hit[2]);
612 track->CalculateClusterWidths(i);
614 Double_t beta = track->GetCrossingAngle(i);
615 track->SetCrossingAngleLUT(i,beta);
617 //if(hit[1]<0 || hit[2]>445)
618 //if(hit[2]<0 || hit[2]>445)
619 //cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
620 //cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
625 cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
628 void AliL3Modeller::CheckForOverlaps(Float_t dangle,Int_t *rowrange)
630 //Flag the tracks that overlap
633 cout<<"Checking for overlaps on "<<fTracks->GetNTracks()<<endl;
636 for(Int_t k=AliL3Transform::GetFirstRow(fPatch); k<=AliL3Transform::GetLastRow(fPatch); k++)
640 if(k < rowrange[0]) continue;
641 if(k > rowrange[1]) break;
643 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
645 AliL3ModelTrack *track1 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
646 if(!track1) continue;
647 if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0) continue;
649 for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
651 AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
652 if(!track2) continue;
653 if(track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0) continue;
655 if(abs((Int_t)rint(track1->GetPadHit(k))-(Int_t)rint(track2->GetPadHit(k))) <= fPadOverlap &&
656 abs((Int_t)rint(track1->GetTimeHit(k))-(Int_t)rint(track2->GetTimeHit(k))) <= fTimeOverlap)
658 if(dangle>0 && fabs(track1->GetCrossingAngleLUT(k) - track2->GetCrossingAngleLUT(k)) < dangle)
661 //cout<<"row "<<k<<" "<<i<<" "<<j<<" "<<track1->GetPadHit(k)<<" "<<track2->GetPadHit(k)<<" "<<fabs(track1->GetCrossingAngleLUT(k) - track2->GetCrossingAngleLUT(k))<<endl;
664 track1->SetOverlap(k,j);
672 cout<<"and there are "<<fTracks->GetNTracks()<<" track left"<<endl;
673 //cout<<"found "<<counter<<" done"<<endl;
677 void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
679 // calculates cluster's width
682 padw = AliL3Transform::GetPadPitchWidth(fPatch);
684 Float_t charge = (Float_t)cl->fCharge;
685 Float_t pad = (Float_t)cl->fPad/charge;
686 Float_t time = (Float_t)cl->fTime/charge;
687 Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
689 //Save the sigmas in pad and time:
691 sigmaY2 = (s2);// + 1./12);//*padw*padw;
693 /*Constants added by offline
696 sigmaY2 = sigmaY2*0.108;
698 sigmaY2 = sigmaY2*2.07;
702 s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
703 timew = AliL3Transform::GetZWidth();
704 sigmaZ2 = (s2);// +1./12);//*timew*timew;
709 Constants added by offline
712 sigmaZ2 = sigmaZ2*0.169;
714 sigmaZ2 = sigmaZ2*1.77;
720 void AliL3Modeller::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
723 AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fRowData;
725 trackID[0]=trackID[1]=trackID[2]=-2;
727 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
729 if(rowPt->fRow < (UInt_t)fCurrentPadRow)
731 AliL3MemHandler::UpdateRowPointer(rowPt);
734 AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
735 for(UInt_t j=0; j<rowPt->fNDigit; j++)
737 Int_t cpad = digPt[j].fPad;
738 Int_t ctime = digPt[j].fTime;
739 if(cpad != pad) continue;
740 if(ctime != time) continue;
741 //if(cpad != pad && ctime != ctime) continue;
742 //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
743 trackID[0] = digPt[j].fTrackID[0];
744 trackID[1] = digPt[j].fTrackID[1];
745 trackID[2] = digPt[j].fTrackID[2];
747 //cout<<"Reading trackID "<<trackID[0]<<endl;
752 void AliL3Modeller::GetTrackID(Int_t /*pad*/,Int_t /*time*/,Int_t */*trackID*/)
754 // Does nothing if do_mc undefined