TPC/TestTPCTrackHits.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.3  2000/11/02 10:22:50  kowal2
  19 Logs added
  20
  21 */
  22
  23 /*
  24   Author : MI
  25   macro to compare TClonesArray hits with interpolated hits
  26   ConvertHits1 read
  27 */
  28
  29 #include "alles.h"
  30 #include "AliObjectArray.h"
  31 #include "AliTPCTrackHits.h"
  32 #include "AliArrayBranch.h"
  33 #include "TestTPCTrackHits.h"
  34
  35 ClassImp(AliTPChitD)
  36
  37 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits,  Bool_t debug, TClonesArray *arrd=0);
  38 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits);
  39
  40 void ConvertHits(const char * benchmark, Bool_t debug)
  41 {
  42   //
  43   //convert - not parametrised hits stored in Clones array
  44   //to
  45
  46
  47   TFile f("galice.root","update");
  48   TClonesArray *arr = new TClonesArray("AliTPChit",100);
  49   TTree * treeCl =(TTree*)f.Get("TreeH0");
  50   TBranch *branch = treeCl->GetBranch("TPC");
  51   AliTPCTrackHits *  myhits = new AliTPCTrackHits;
  52   branch->SetAddress(&arr);
  53   //particle
  54   TTree * treeP = (TTree*)f.Get("TreeK0");
  55   TBranch *branchP = treeP->GetBranch("Particles");
  56   TClonesArray *arrp = new TClonesArray("TParticle",100);
  57   branchP->SetAddress(&arrp);
  58   branchP->GetEvent(0);
  59   //
  60   TFile f2("treeh.root","recreate");
  61   f2.SetCompressionLevel(10);
  62   TTree * treeh2 = new TTree("TreeTPCH","TreeTPCH");
  63   //treeh2->AliBranch("TPC","AliTPCTrackHits", &myhits,4096);
  64   treeh2->GetListOfBranches()->Add(new AliObjectBranch("TPC","AliTPCTrackHits",
  65          &myhits,treeh2,4096,1));
  66
  67   TFile f3("treehcl.root","recreate");
  68   f3.SetCompressionLevel(2);
  69   TTree * treeh3 = new TTree("TreeTPCH","TreeTPCH");
  70   treeh3->Branch("TPC", &arr,64000,100);
  71
  72   gBenchmark->Start(benchmark);
  73   Int_t trsize = (Int_t)treeCl->GetEntries();
  74   for (Int_t i=0;i<300;i++){
  75     arr->Clear();
  76     Int_t size = branch->GetEvent(i);
  77     //if (size>0){
  78     if ((size>0)&& arr->GetEntriesFast()>0) {
  79       f.cd();
  80       myhits->Clear();
  81       myhits =MakeTrack(arr,arrp,myhits);
  82       CompareHits(arr,myhits,debug);
  83       f2.cd();
  84       treeh2->Fill();
  85       f3.cd();
  86       treeh3->Fill();
  87       if ((i%10)==0) cout<<i<<"\n";
  88     }
  89   }
  90   f3.cd();
  91   treeh3->Write();
  92   f2.cd();
  93   treeh2->Write();
  94   f2.Close();
  95   f.Close();
  96   delete myhits;
  97   gBenchmark->Stop(benchmark);
  98   gBenchmark->Show(benchmark);
  99 }
 100
 101
 102
 103 void CompareHits(const char * benchmark, Bool_t debug)
 104 {
 105   //compare persistent hits
 106   TFile f2("treeh.root");
 107   TTree * treeh2 = (TTree*)f2.Get("TreeTPCH");
 108   AliTPCTrackHits *  myhits = new AliTPCTrackHits ;
 109   AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC");
 110   mybranch->SetAddress(&myhits);
 111
 112
 113   TClonesArray *arr = new TClonesArray("AliTPChit",100);
 114   TFile f3("treehcl.root");
 115   TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
 116   TBranch *branch = treeh3->GetBranch("TPC");
 117   branch->SetAddress(&arr);
 118
 119   cout<<"Lets go!\n";
 120   gBenchmark->Start(benchmark);
 121   Int_t trsize = (Int_t)treeh3->GetEntries();
 122   for (Int_t i=0;i<300;i++){
 123     Int_t size = branch->GetEvent(i);
 124     mybranch->GetEvent(i);
 125     //if (arr){
 126     if ((arr->GetEntriesFast()>0) && size>0) {
 127       CompareHits(arr,myhits,debug);
 128       if ((i%10)==0) cout<<i<<"\n";
 129     }
 130   }
 131   gBenchmark->Stop(benchmark);
 132   gBenchmark->Show(benchmark);
 133 }
 134
 135
 136 void CompareHitsG(const char * benchmark, Bool_t debug)
 137 {
 138   //compare persistent hits
 139   TFile f2("galice.root");
 140   TTree * treeh2 = (TTree*)f2.Get("TreeH0");
 141   AliTPCTrackHits *  myhits = new AliTPCTrackHits ;
 142   AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC2");
 143   mybranch->SetAddress(&myhits);
 144
 145
 146   TClonesArray *arr = new TClonesArray("AliTPChit",100);
 147   //TFile f3("treehcl.root");
 148   //TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
 149   TBranch *branch = treeh2->GetBranch("TPC");
 150   branch->SetAddress(&arr);
 151
 152   TFile f3("treehdelta.root","recreate");
 153   f3.SetCompressionLevel(2);
 154   TTree * treeh3 = new TTree("DelataH","DeltaH");
 155   TClonesArray *arrd = new TClonesArray("AliTPChitD",100);
 156   treeh3->Branch("TPC", &arrd,64000,100);
 157
 158   cout<<"Lets go!\n";
 159   gBenchmark->Start(benchmark);
 160   Int_t trsize = treeh2->GetEntries();
 161   for (Int_t i=0;i<trsize;i++){
 162     Int_t size = branch->GetEvent(i);
 163     mybranch->GetEvent(i);
 164     //if (arr){
 165     if ((arr->GetEntriesFast()>0) && size>0) {
 166       CompareHits(arr,myhits,debug,arrd);
 167     }
 168     if ((i%10)==0) cout<<i<<"\n";
 169     treeh3->Fill();
 170   }
 171   treeh3->Write();
 172   f3.Close();
 173   gBenchmark->Stop(benchmark);
 174   gBenchmark->Show(benchmark);
 175 }
 176
 177
 178
 179
 180 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits)
 181 {
 182   //
 183   //make track wit hits  according
 184   AliTPChit * hit;
 185   //  AliTPCTrackHits * myhits= new AliTPCTrackHits;
 186   myhits->SetHitPrecision(0.002);
 187   myhits->SetStepPrecision(0.003);
 188   myhits->SetMaxDistance(100);
 189   myhits->FlushHitStack(kTRUE);
 190   for (Int_t i=0;i<arr->GetEntriesFast();i++){
 191     hit = (AliTPChit*)arr->At(i);
 192     if (hit){
 193       TParticle *p = (TParticle*)arrp->At(hit->GetTrack());
 194       Float_t momentum = TMath::Sqrt(p->Px()*p->Px()+p->Py()*p->Py());
 195       Float_t ran= 100.*gRandom->Rndm();
 196       if (ran<1.)  myhits->FlushHitStack(kTRUE);
 197       if (momentum<0.01) {//if small momentum - not precise
 198         myhits->SetHitPrecision(0.05);
 199         myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
 200                              hit->Z(), hit->fQ+1000);
 201       }
 202       else {
 203         myhits->SetHitPrecision(0.002);
 204         myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
 205                              hit->Z(), hit->fQ);
 206       }
 207     }
 208   }
 209   myhits->FlushHitStack();
 210   return myhits;
 211 }
 212
 213
 214
 215 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd)
 216 {
 217   //
 218   // if debug option  kTRUE
 219   // compare hits and write result  to the stdoutput
 220   AliTPChit * hit, *hit2;
 221   if (arrd) arrd->Clear();
 222
 223   for (Int_t i=0;i<arr->GetEntriesFast();i++){
 224     hit = (AliTPChit*)arr->At(i);
 225     if (hit){
 226       if (i==0) myhits->First();
 227       else myhits->Next();
 228       hit2 = myhits->GetHit();
 229       if (!hit2) {
 230         hit2=0;
 231         if (hit) cout<<"Error _ hits "<<i<<"didn't find\n";
 232       }
 233
 234       if (hit&&hit2){
 235         // if (hit){
 236
 237         AliTrackHitsParam *param= myhits->GetParam();
 238         AliHitInfo * info = myhits->GetHitInfo();
 239         //
 240         if (arrd) {
 241           TClonesArray &larrd = *arrd;
 242           AliTPChitD * h = new(larrd[i]) AliTPChitD;
 243           h->SetX(hit->X());
 244           h->SetY(hit->Y());
 245           h->SetZ(hit->Z());
 246           h->SetTrack(hit->GetTrack());
 247           h->fQ = hit->fQ;
 248           h->fSector = hit->fSector;
 249           AliTPChit * hitd = h->GetDelta();
 250           hitd->SetX(hit->X()-hit2->X());
 251           hitd->SetY(hit->Y()-hit2->Y());
 252           hitd->SetZ(hit->Z()-hit2->Z());
 253           hitd->SetTrack(hit->GetTrack()-hit2->GetTrack());
 254           hitd->fQ = hit->fQ-hit2->fQ;
 255           hitd->fSector = hit->fSector-hit2->fSector;
 256         }
 257
 258         if (debug){
 259           Float_t dd =
 260             TMath::Sqrt(
 261                         (hit->X()-hit2->X())*(hit->X()-hit2->X())+
 262                         (hit->Y()-hit2->Y())*(hit->Y()-hit2->Y())+
 263                         (hit->Z()-hit2->Z())*(hit->Z()-hit2->Z()));
 264           printf("C1\t%d\t%d\t%d\t%d\t",
 265                  hit->fSector,hit2->fSector,hit->GetTrack(),hit2->GetTrack());
 266           printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
 267                  hit->X(),hit2->X(), hit->Y(),hit2->Y(),hit->Z(),hit2->Z());
 268           printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
 269                  dd, param->fR,param->fZ,param->fFi,param->fAn,
 270                  param->fAd,param->fTheta);
 271           printf("%d\t%d\t%d\n",
 272                  (Int_t)info->fHitDistance, (Int_t)hit->fQ, (Int_t)hit2->fQ);
 273         }
 274       }
 275     }
 276   }
 277 }
 278
 279
 280
 281
 282 //Filter for paw visualisation of results
 283 //
 284 //sed filter
 285  //sed '/*/d' dd.txt | sed '/C/d'|  cat -n > out.txt
 286  // sed -n '/C2/p' dd.txt | sed 's/C1//' | cat -n >out2.txt
 287  // sed -n '/C3/p' dd.txt | sed 's/C2//' | cat -n >out3.txt
 288
 289 //filter
 290 //myfilter C1 dd.txt >out1.txt
 291 //myfilter C2 dd.txt >out2.txt
 292 //myfilter C3 dd.txt >out3.txt
 293 // sed -n 1,50000p
 294 // sed -n 1,50000p dd.txt | myfilter C1  | sed -n 1,50000p >out1.txt
 295
 296 //myfilter C1 dd.txt | sed -n 1,50000p >out1.txt
 297 //myfilter C2 dd.txt | sed -n 1,50000p >out2.txt
 298 //myfilter C3 dd.txt | sed -n 1,50000p >out3.txt
 299 /*
 300 paw visualisation
 301  Nt/Create 1 'Test of Match' 21 ! ! event v1 v2 t1 t2  x1 x2 y1  y2 z1 z2 dd r z fi an ad theta dist q1 q2
 302   Nt/Read 1 out1.txt
 303  nt
 304  Nt/Create 2 'Test of Match2' 13 ! ! event stacksize i r  z fi alfa theta dr1 ddr ddz ddrfi dd
 305   Nt/Read 2 out2.txt
 306
 307  Nt/Create 3 'Test of Match2' 2 ! ! event stacksize
 308   Nt/Read 3 out3.txt
 309 */
 310
 311 void   Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
 312             Double_t fSumX,  Double_t fSumX2, Double_t fSumX3,
 313             Double_t fSumX4, Int_t n,
 314             Double_t &a, Double_t &b, Double_t &c)
 315 {
 316   //
 317   //recalc parameters not fixing origin point
 318   Double_t det =
 319     n* (fSumX2*fSumX4-fSumX3*fSumX3) -
 320     fSumX*      (fSumX*fSumX4-fSumX3*fSumX2)+
 321     fSumX2*     (fSumX*fSumX3-fSumX2*fSumX2);
 322
 323   if (TMath::Abs(det)>0){
 324     a =
 325       (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
 326        fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
 327        fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det;
 328     b=
 329       (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
 330       fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
 331       fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
 332     c=
 333       (n*(fSumX2*fSumYX2-fSumYX*fSumX3)-
 334        fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
 335        fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;
 336     cout<<a<<"\t"<<b<<"\t"<<c<<"\n";
 337   }
 338 }
 339
 340 void TestFit(Float_t a, Float_t b, Float_t c, Int_t n)
 341 {
 342   Double_t fSumY,fSumYX,fSumYX2,fSumX, fSumX2,fSumX3, fSumX4;
 343   fSumY = fSumYX = fSumYX2 = fSumX =  fSumX2 = fSumX3 =  fSumX4 =0;
 344   Double_t a2,b2,c2;
 345   for (Int_t i=0;i<n; i++){
 346     Float_t  x = gRandom->Rndm();
 347     Float_t  y = a+b*x+c*x*x;
 348     fSumY+=y;
 349     fSumYX+=y*x;
 350     fSumYX2+=y*x*x;
 351     fSumX += x;
 352     fSumX2 += x*x;
 353     fSumX3 += x*x*x;
 354     fSumX4 += x*x*x*x;
 355   }
 356   Fit2(fSumY,fSumYX,fSumYX2, fSumX,fSumX2,fSumX3,fSumX4,n,a2,b2,c2);
 357 }