]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/TPCDigits2Clusters.C
Class for fixed point operations.
[u/mrichter/AliRoot.git] / TPC / TPCDigits2Clusters.C
index ed3809e985e6fd2e4777b85e018c1e8df5e83219..580cd6c85c60200f553a260363404551fdf44cfb 100644 (file)
-void TPCDigits2Clusters() {
-// Dynamically link some shared libs
-   if (gClassTable->GetID("AliRun") < 0) {
-      gSystem->Load("libGeant3Dummy.so");      // a dummy version of Geant3
-      gSystem->Load("PHOS/libPHOSdummy.so");   // the standard Alice classes 
-      gSystem->Load("libgalice.so");           // the standard Alice classes 
-   }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
-   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
-   if (file) file->Close();
-   if (!file) file = new TFile("galice.root");
-
-// Get AliRun object from file or create it if not on file
-   if (!gAlice) {
-      gAlice = (AliRun*)file->Get("gAlice");
-      if (gAlice) printf("AliRun object found on file\n");
-      if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
-   }
-
-   gAlice->GetEvent(0);
-
-   AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
-   TPC->Digits2Clusters();
-   TClonesArray *c=TPC->Clusters();
-   int n=c->GetEntriesFast();
-   cout<<"Number of clusters "<<n<<endl;
-
-   TPolyMarker3D *pm=new TPolyMarker3D(n);
-   for (int i=0; i<n; i++) {
-       AliTPCcluster *cl=(AliTPCcluster *)c->UncheckedAt(i);
-       pm->SetPoint(i,cl->fX,cl->fY,cl->fZ);
-   }
-
-   c1=new TCanvas("c1", "Cluster display",0,0,575,750);
-   TView *v=new TView(1);
-   v->SetRange(-430,-560,-430,430,560,1710);
-
-   c1->Clear();
-   c1->SetFillColor(1);
-   c1->SetTheta(90.);
-   c1->SetPhi(0.);
-
-   pm->SetMarkerSize(1);
-   pm->SetMarkerColor(2);
-   pm->SetMarkerStyle(1);
-   pm->Draw();
-
-   gAlice->GetGeometry()->Draw("same");
+
+/*
+#include "alles.h"
+extern AliTPCParam *gTPCParam;
+extern AliTPCClustersArray * gCalcClusters;
+extern AliTPCClustersArray * gDifClusters;
+AliTPCDigitsArray * GetDigitsArray();
+  AliTPCClustersArray * GetExactClustersArray();
+AliTPCClustersArray *  GetCalcClustersArray(Bool_t newtree=kFALSE);
+AliTPCClustersArray *  GetDifClustersArray(Bool_t newtree=kFALSE);
+TClonesArray * CompareClusters(TClonesArray *calcclusters,TClonesArray *exactclusters, 
+                               TClonesArray *diffclusters, Float_t shx, Float_t shy);
+*/
+
+void   TPCDigits2Clusters(AliTPCDigitsArray * digarr, AliTPCClustersArray *calcarr,
+                         AliTPCClustersArray *exarr=0, AliTPCClustersArray *difarr=0)
+{
+  //
+  //calculate clusters position from digits information
+  //
+
+  
+  //initialisation of cluster finder
+  AliClusterFinder cf;
+  cf.SetThreshold(gTPCParam->GetZeroSup());
+  cf.SetDetectorParam(gTPCParam);
+  cf.SetBFit(kFALSE);
+  cf.GetMinuit()->SetPrintLevel(-1);
+  cf.GetMinuit()->SetMaxIterations(20);
+  
+  if (digarr ==0) digarr = GetDigitsArray();
+  if ( (digarr == 0) || (digarr->GetTree()==0)) return;  
+  if (exarr==0) exarr =GetExactClustersArray();
+  if ((exarr!=0)&&(difarr==0)) difarr= GetDifClustersArray();
+
+  if (calcarr==0) calcarr = GetCalcClustersArray(kTRUE);
+  //loop over all writen digits segments
+  Int_t nsegment = (Int_t)digarr->GetTree()->GetEntries();
+  Int_t segment;
+  for (segment = 0; segment<nsegment;segment++){
+    //load segment with index (TTree internal)  number segment
+    AliSimDigits *digrow= (AliSimDigits*)digarr->LoadEntry(segment);
+    Int_t sector,row;
+    ((AliTPCParam*)digarr->GetParam())->AdjustSectorRow(digrow->GetID(),sector,row);
+    
+    AliTPCClustersRow * clrow = calcarr->CreateRow(sector,row);
+    AliTPCClustersRow * difrow= difarr->CreateRow(sector,row);    
+    AliTPCClustersRow * exrow =  exarr->GetRow(sector,row);
+    if (exrow==0) exrow = exarr->LoadRow(sector,row);
+    TH2F *his = (TH2F*)digrow->GenerHisto();
+    digrow->ExpandTrackBuffer();
+    if (his==0) return;
+
+    //set current index for cluster finder
+    Int_t  index[3]= {0,sector,row};
+    cf.GetHisto(his);
+    cf.SetDetectorIndex(index);    
+    cf.FindPeaks3(clrow->GetArray());
+
+    //get cluster tracks ID
+    Int_t ncl = clrow->GetArray()->GetEntriesFast();
+    for (Int_t icl = 0 ;icl<ncl; icl++)
+      {
+       AliDigitCluster *  cl =(AliDigitCluster*)clrow->GetArray()->UncheckedAt(icl);
+       Int_t i = (Int_t)cl->fMaxX;
+       Int_t j = (Int_t)cl->fMaxY;             
+       cl->fTracks[0]=digrow->GetTrackIDFast(i,j,0);
+       cl->fTracks[1]=digrow->GetTrackIDFast(i,j,1);
+       cl->fTracks[2]=digrow->GetTrackIDFast(i,j,2);   
+      }
+    Float_t shx = gTPCParam->GetZOffset()/gTPCParam->GetZWidth()+0.5;
+    Float_t shy= 0.5;
+    CompareClusters(clrow->GetArray(),exrow->GetArray(),difrow->GetArray(),shx,shy);        
+    calcarr->StoreRow(sector,row);
+    difarr->StoreRow(sector,row);
+    digarr->ClearRow(sector,row);
+  }    
+  //STORING RESULTS
+ char treeName[100];
+ sprintf(treeName,"TreeCCalc_%s",gTPCParam->GetTitle());
+ calcarr->GetTree()->Write(treeName);
+ if (difarr!=0){
+   char treeName[100];
+   sprintf(treeName,"TreeCDif_%s",gTPCParam->GetTitle());
+   difarr->GetTree()->Write(treeName);
+ }
+
 }
+
+
+AliTPCClustersArray *  GetCalcClustersArray(Bool_t newtree=kFALSE, const char* name=0) 
+{
+  //
+  //construct AliTPCClusters Array object
+  //
+  AliTPCClustersArray * arr=0;
+  //  if ( (gAlice!=0) && (gAlice->GetDetector("TPC")!=0) ) {    
+  //  arr = ((AliTPC*)gAlice->GetDetector("TPC"))->GetClustersArray();
+  //  if (arr!=0) arr->Update();
+  //}
+  if (arr==0) {    
+    arr = new AliTPCClustersArray;
+    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+    if (file==0) file = TFile::Open("rfio:galice.root","update");   
+    arr->SetClusterType("AliDigitCluster");         
+    arr->Setup(gTPCParam);
+    cout<<"Update status : "<<arr->Update()<<"\n";
+    char treeName[100];
+    if (name==0)
+      sprintf(treeName,"TreeCExact_%s",gTPCParam->GetTitle());
+    else
+      sprintf(treeName,"TreeCCalc_%s",gTPCParam->GetTitle());    
+    if (newtree!=kTRUE){
+      cout<<"Connect tree status : "<<arr->ConnectTree(treeName)<<"\n";
+    }
+    else {
+      arr->MakeTree();
+    }
+  }
+  gCalcClusters = arr;
+  return arr;
+}
+
+
+AliTPCClustersArray *  GetDifClustersArray(Bool_t newtree=kFALSE, const char* name=0) 
+{
+  //
+  //construct AliTPCClusters Array object
+  //
+  AliTPCClustersArray * arr=0;
+  //  if ( (gAlice!=0) && (gAlice->GetDetector("TPC")!=0) ) {    
+  //  arr = ((AliTPC*)gAlice->GetDetector("TPC"))->GetClustersArray();
+  //  if (arr!=0) arr->Update();
+  //}
+  if (arr==0) {    
+    arr = new AliTPCClustersArray;
+    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+    if (file==0) file = TFile::Open("rfio:galice.root","update");   
+    arr->SetClusterType("AliDifCluster");         
+    arr->Setup(gTPCParam);
+    cout<<"Update status : "<<arr->Update()<<"\n";
+    char treeName[100];
+    if (name==0) 
+      sprintf(treeName,"TreeCExact_%s",gTPCParam->GetTitle());
+    else
+      sprintf(treeName,"TreeCDif_%s",gTPCParam->GetTitle());
+    if (newtree!=kTRUE){
+      cout<<"Connect tree status : "<<arr->ConnectTree(treeName)<<"\n";
+    }
+    else {
+      arr->MakeTree();
+    }
+  }
+  gDifClusters  = arr; 
+  return arr;
+}
+
+
+
+TClonesArray * CompareClusters(TClonesArray *calcclusters,TClonesArray *exactclusters, 
+                               TClonesArray *diffclusters, Float_t shx, Float_t shy)
+{
+  if (calcclusters==0) return 0;
+  if (exactclusters==0) return 0;
+  if (diffclusters==0) return 0;
+  TClonesArray & diff=*diffclusters;
+   
+   
+  Int_t nclusters2 = calcclusters->GetEntriesFast();
+  Int_t nclusters1 = exactclusters->GetEntriesFast();
+  
+  for(Int_t i=0; i<nclusters1; i++){
+    AliCluster * cl1 = (AliCluster*)exactclusters->UncheckedAt(i);
+    AliDigitCluster *cl2;
+    AliDifCluster diffc;
+    Int_t index=-1;
+    Float_t dx=10000;
+    Float_t dy=10000;
+    for (Int_t j=0; j<nclusters2;j++){
+      cl2 = (AliDigitCluster*)calcclusters->UncheckedAt(j);
+      if ( (cl2!=0)&& (cl1!=0)){ 
+       Float_t ddx = (cl2->fX-shx)-cl1->fX;
+       Float_t ddy = (cl2->fY-shy)-cl1->fY;
+       if ((ddx*ddx+ddy*ddy)<(dx*dx+dy*dy)){
+         dx=ddx;
+         dy=ddy;
+         index=j;
+       }      
+      }
+    }
+    
+    cl2 = (AliDigitCluster*)calcclusters->UncheckedAt(index);
+    if (cl2!=0){
+      diffc.fDx      =dx;
+      diffc.fDy      =dy;
+      
+      diffc.fAngleX  = cl1->fSigmaX2;
+      diffc.fAngleY  = cl1->fSigmaY2; 
+      diffc.fTracks[0] = cl2->fTracks[0];
+      diffc.fTracks[1] = cl2->fTracks[1];
+      diffc.fTracks[2] = cl2->fTracks[2];
+      diffc.fGenerTrack= cl1->fTracks[0];
+
+      diffc.fX       = cl2->fX;
+      diffc.fY       = cl2->fY;
+      diffc.fNx      = cl2->fNx;
+      diffc.fNy      = cl2->fNy;
+      diffc.fQ       = cl2->fQ;    
+      diffc.fOrigQ   = cl1->fQ;
+      diffc.fSigmaX2 = cl2->fSigmaX2;
+      diffc.fSigmaY2 = cl2->fSigmaY2;
+      diffc.fArea    = cl2->fArea;
+      diffc.fMax     = cl2->fMax;
+    }
+    else cout<<"pici[ici/n";
+    
+    new(diff[i]) AliDifCluster(diffc);
+  }
+  return diffclusters;
+}