HLT/hough/test.C

   1 void test(char *file="/prog/alice/data/Rawdata/1_patch/pp/event_0/",bool bin=true)
   2 {
   3   AliL3Logger l;
   4 //  l.UnSet(AliL3Logger::kDebug);
   5 //  l.UnSet(AliL3Logger::kAll);
   6   l.Set(AliL3Logger::kAll);
   7   //l.UseStdout();
   8   l.UseStream();
   9
  10   Int_t fNEtaSegments = 1;
  11
  12   Char_t histname[50];
  13   Int_t i;
  14
  15   int slice=1,patch=0;
  16
  17   fHoughTransformer = new AliL3HoughTransformer(slice,patch,fNEtaSegments);
  18   fMemHandler = new AliL3FileHandler();
  19   fMaxFinder = new AliL3HoughMaxFinder("KappaPhi");
  20
  21   fTransform = new AliL3Transform();
  22
  23   fHoughTransformer->CreateHistograms(64,-0.012,0.012,64,-0.35,0.35);
  24   fHoughTransformer->SetThreshold(10);
  25   fTracks = new AliL3TrackArray("AliL3HoughTrack");
  26
  27   TH2F *road = new TH2F("road","",250,0,250,250,-125,125);
  28   TH2F *peaks = new TH2F("peaks","",64,-0.006,0.006,64,-0.26,0.26);
  29   peaks->SetMarkerStyle(3);
  30   peaks->SetMarkerColor(2);
  31   road->SetMarkerStyle(5);
  32   road->SetMarkerColor(2);
  33
  34   Char_t name[256];
  35
  36   UInt_t ndigits=0;
  37   AliL3DigitRowData *digits =0;
  38
  39   if(bin) //read data from binary files
  40     {
  41       Char_t fPath[256];
  42       strcpy(fPath,file);
  43       fMemHandler->Free();
  44       sprintf(name,"%sdigits_%d_%d.raw",fPath,slice,patch);
  45       fMemHandler->SetBinaryInput(name);
  46       digits = (AliL3DigitRowData *)fMemHandler->CompBinary2Memory(ndigits);
  47       fMemHandler->CloseBinaryInput();
  48     }
  49   else //read data from root file
  50     {
  51       const Int_t NRows[6][2] = {{0,31},{32,63},{64,91},{92,119},{120,143},{144,175}};
  52       const int nr[2] = {0,175};
  53       fMemHandler->Free();
  54       fMemHandler->SetAliInput(file);
  55       fMemHandler->Init(slice,patch,NRows[patch]);
  56       fMemHandler->Init(fTransform);
  57       digits=(AliL3DigitRowData *)fMemHandler->AliDigits2Memory(ndigits);
  58       fMemHandler->CloseAliInput();
  59     }
  60
  61   fHoughTransformer->SetInputData(ndigits,digits);
  62   fEval = new AliL3HoughEval();
  63   fEval->InitTransformer(fHoughTransformer);
  64   //fEval->SetNumOfRowsToMiss(0);
  65   //fEval->SetNumOfPadsToLook(2);
  66   //fEval->RemoveFoundTracks();
  67   AliL3HoughTrack *track;
  68   AliL3Histogram *histo;
  69   Int_t good_count;
  70   int eind=0;
  71   while(1)
  72     {
  73       printf("Transforming\n");
  74       double init = AliL3Benchmark::GetCpuTime();
  75       fHoughTransformer->TransformCircle();
  76       double final = AliL3Benchmark::GetCpuTime();
  77       printf("done in %f ms\n",(final-init)*1000);
  78       good_count=0;
  79       for(Int_t e=0; e<fNEtaSegments; e++)
  80         {
  81           histo = (AliL3Histogram*)fHoughTransformer->GetHistogram(e);
  82           if(!histo) continue;
  83
  84           fMaxFinder->SetHistogram(histo);
  85           Int_t n=10;
  86           Int_t weight[10];
  87           Float_t x[10];
  88           Float_t y[10];
  89           fMaxFinder->FindPeak1(x,y,weight,n);
  90           track = new AliL3HoughTrack();
  91           track->SetTrackParameters(x[0],y[0],1);
  92
  93           //      if(!fEval->LookInsideRoad(track,e))
  94           // continue;
  95
  96
  97           if(e==eind)
  98             {
  99               for(int t=0; t<176; t++)
 100                 {
 101                   if(t%10) continue;
 102                   float xyz_tr[3];
 103                   track->GetCrossingPoint(t,xyz_tr);
 104                   road->Fill(xyz_tr[0],xyz_tr[1],1);
 105                 }
 106             }
 107           delete track;
 108
 109         }
 110
 111       break;
 112       if(good_count==0)
 113         break;
 114     }
 115
 116
 117   image = new AliL3Histogram("image","",250,0,250,250,-125,125);
 118   fEval->DisplayEtaSlice(eind,image);
 119
 120   c1 = new TCanvas("c1","",1000,500);
 121   c1->Divide(2);
 122   c1->cd(1);
 123   histo = (AliL3Histogram*)fHoughTransformer->GetHistogram(eind);
 124   if(!histo)
 125     {printf("No histogram\n"); return;}
 126   histo->Draw("box");
 127
 128   peaks->Draw("same");
 129   c1->cd(2);
 130   image->Draw("");
 131   road->Draw("same");
 132 }
 133
 134 void process(char *digitfile="/prog/alice/data/Rawdata/6_patch/hg_42105_s1-3/",bool bin=true,char *trackfile="good_tracks")
 135 {
 136   AliL3Logger l;
 137   //  l.UnSet(AliL3Logger::kDebug);
 138   //  l.UnSet(AliL3Logger::kAll);
 139   //l.Set(AliL3Logger::kError);
 140   //l.UseStdout();
 141   l.UseStream();
 142
 143   transform = new AliL3Transform();
 144
 145   double torad = 3.1415/180;
 146   int slice=1;
 147   a = new AliL3Hough(digitfile,bin,100);
 148
 149   a->ReadData(slice);
 150   a->Transform();
 151   a->AddAllHistograms();
 152   a->FindTrackCandidates();
 153   a->Evaluate(2);
 154   a->WriteTracks("../comp/");
 155
 156   //a->MergePatches();
 157
 158   //a->MergeInternally();
 159
 160   //tracks = (AliL3TrackArray*)(a->GetInterMerger())->GetOutTracks();
 161   //tracks = (AliL3TrackArray*)(a->GetMerger())->GetOutTracks();
 162   tracks = (AliL3TrackArray*)a->GetTracks(0);
 163   //a->GetEval(0)->CompareMC(tracks,"good_tracks_hg4000_s1");
 164
 165   float xyz[3];
 166   int row=0;
 167   for(int i=0; i<tracks->GetNTracks(); i++)
 168     {
 169       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
 170       if(!track) continue;
 171       track->GetCrossingPoint(row,xyz);
 172       transform->Local2Raw(xyz,1,row);
 173       cout<<"Recon. pad "<<(int)xyz[1]<<" time "<<(int)xyz[2]<<endl;
 174       cout<<"Pt "<<track->GetPt()<<" phi "<<track->GetPhi0()<<" charge "<<track->GetCharge()<<endl;
 175     }
 176
 177   delete transform;
 178   return;
 179
 180 }
 181
 182 void cf()
 183 {
 184     AliL3Logger l;
 185   //  l.UnSet(AliL3Logger::kDebug);
 186   //  l.UnSet(AliL3Logger::kAll);
 187   //l.Set(AliL3Logger::kError);
 188   //l.UseStdout();
 189   l.UseStream();
 190
 191   AliL3FileHandler *file = new AliL3FileHandler();
 192   file->SetBinaryInput("tracks.raw");
 193   fTracks = new AliL3TrackArray("AliL3HoughTrack");
 194   file->Binary2TrackArray(fTracks);
 195   cout<<"Ntracks "<<fTracks->GetNTracks()<<endl;
 196   file->CloseBinaryInput();
 197   delete file;
 198
 199
 200   return;
 201   /*
 202   a = new AliL3ClusterFinder();
 203   a->LoadData();
 204   a->Process();
 205   */
 206 }
 207
 208
 209 void HistogramParticles(char *trackfile)
 210 {
 211   struct GoodTrack
 212   {
 213
 214     Int_t label;
 215     Double_t eta;
 216     Int_t code;
 217     Double_t px,py,pz;
 218     Double_t pt;
 219     Int_t nhits;
 220   };
 221   struct GoodTrack goodtracks[15000];
 222   Int_t nt=0;
 223   ifstream in(trackfile);
 224   if(in)
 225     {
 226       printf("Reading good tracks from file %s\n",trackfile);
 227       while (in>>goodtracks[nt].label>>goodtracks[nt].code>>
 228              goodtracks[nt].px>>goodtracks[nt].py>>goodtracks[nt].pz>>
 229              goodtracks[nt].pt>>goodtracks[nt].eta>>goodtracks[nt].nhits)
 230         {
 231           nt++;
 232           if (nt==15000)
 233             {
 234               cerr<<"Too many good tracks"<<endl;
 235               break;
 236             }
 237         }
 238     }
 239
 240   int fNEtaSegments = 150;
 241   Double_t etaslice = (0.9 - 0)/fNEtaSegments;
 242   Int_t *particles = new Int_t[fNEtaSegments];
 243   for(Int_t i=0; i<fNEtaSegments; i++)
 244     particles[i]=0;
 245
 246   for(Int_t i=0; i<nt; i++)
 247     {
 248       //if(goodtracks[i].nhits < 150) continue;
 249       if(goodtracks[i].pt < 0.5) continue;
 250       Int_t particleindex = (Int_t)(goodtracks[i].eta/etaslice);
 251       if(particleindex < 0 || particleindex >= fNEtaSegments) continue;
 252       particles[particleindex]++;
 253     }
 254
 255   hist = new TH1F("hist","",21,0,20);
 256   for(int i=0; i<fNEtaSegments; i++)
 257     hist->Fill(particles[i]);
 258
 259
 260   c1 = new TCanvas("c1","",2);
 261   hist->Draw();
 262
 263 }