TPC/TPCDigits2Clusters.C

   1
   2 /*
   3 #include "alles.h"
   4 extern AliTPCParam *gTPCParam;
   5 extern AliTPCClustersArray * gCalcClusters;
   6 extern AliTPCClustersArray * gDifClusters;
   7 AliTPCDigitsArray * GetDigitsArray();
   8   AliTPCClustersArray * GetExactClustersArray();
   9 AliTPCClustersArray *  GetCalcClustersArray(Bool_t newtree=kFALSE);
  10 AliTPCClustersArray *  GetDifClustersArray(Bool_t newtree=kFALSE);
  11 TClonesArray * CompareClusters(TClonesArray *calcclusters,TClonesArray *exactclusters,
  12                                 TClonesArray *diffclusters, Float_t shx, Float_t shy);
  13 */
  14
  15 void   TPCDigits2Clusters(AliTPCDigitsArray * digarr, AliTPCClustersArray *calcarr,
  16                           AliTPCClustersArray *exarr=0, AliTPCClustersArray *difarr=0)
  17 {
  18   //
  19   //calculate clusters position from digits information
  20   //
  21
  22
  23   //initialisation of cluster finder
  24   AliClusterFinder cf;
  25   cf.SetThreshold(gTPCParam->GetZeroSup());
  26   cf.SetDetectorParam(gTPCParam);
  27   cf.SetBFit(kFALSE);
  28   cf.GetMinuit()->SetPrintLevel(-1);
  29   cf.GetMinuit()->SetMaxIterations(20);
  30
  31   if (digarr ==0) digarr = GetDigitsArray();
  32   if ( (digarr == 0) || (digarr->GetTree()==0)) return;
  33   if (exarr==0) exarr =GetExactClustersArray();
  34   if ((exarr!=0)&&(difarr==0)) difarr= GetDifClustersArray();
  35
  36   if (calcarr==0) calcarr = GetCalcClustersArray(kTRUE);
  37   //loop over all writen digits segments
  38   Int_t nsegment = (Int_t)digarr->GetTree()->GetEntries();
  39   Int_t segment;
  40   for (segment = 0; segment<nsegment;segment++){
  41     //load segment with index (TTree internal)  number segment
  42     AliSimDigits *digrow= (AliSimDigits*)digarr->LoadEntry(segment);
  43     Int_t sector,row;
  44     ((AliTPCParam*)digarr->GetParam())->AdjustSectorRow(digrow->GetID(),sector,row);
  45
  46     AliTPCClustersRow * clrow = calcarr->CreateRow(sector,row);
  47     AliTPCClustersRow * difrow= difarr->CreateRow(sector,row);
  48     AliTPCClustersRow * exrow =  exarr->GetRow(sector,row);
  49     if (exrow==0) exrow = exarr->LoadRow(sector,row);
  50     TH2F *his = (TH2F*)digrow->GenerHisto();
  51     digrow->ExpandTrackBuffer();
  52     if (his==0) return;
  53
  54     //set current index for cluster finder
  55     Int_t  index[3]= {0,sector,row};
  56     cf.GetHisto(his);
  57     cf.SetDetectorIndex(index);
  58     cf.FindPeaks3(clrow->GetArray());
  59
  60     //get cluster tracks ID
  61     Int_t ncl = clrow->GetArray()->GetEntriesFast();
  62     for (Int_t icl = 0 ;icl<ncl; icl++)
  63       {
  64         AliDigitCluster *  cl =(AliDigitCluster*)clrow->GetArray()->UncheckedAt(icl);
  65         Int_t i = (Int_t)cl->fMaxX;
  66         Int_t j = (Int_t)cl->fMaxY;
  67         cl->fTracks[0]=digrow->GetTrackIDFast(i,j,0);
  68         cl->fTracks[1]=digrow->GetTrackIDFast(i,j,1);
  69         cl->fTracks[2]=digrow->GetTrackIDFast(i,j,2);
  70       }
  71     Float_t shx = gTPCParam->GetZOffset()/gTPCParam->GetZWidth()+0.5;
  72     Float_t shy= 0.5;
  73     CompareClusters(clrow->GetArray(),exrow->GetArray(),difrow->GetArray(),shx,shy);
  74     calcarr->StoreRow(sector,row);
  75     difarr->StoreRow(sector,row);
  76     digarr->ClearRow(sector,row);
  77   }
  78   //STORING RESULTS
  79  char treeName[100];
  80  sprintf(treeName,"TreeCCalc_%s",gTPCParam->GetTitle());
  81  calcarr->GetTree()->Write(treeName);
  82  if (difarr!=0){
  83    char treeName[100];
  84    sprintf(treeName,"TreeCDif_%s",gTPCParam->GetTitle());
  85    difarr->GetTree()->Write(treeName);
  86  }
  87
  88 }
  89
  90
  91 AliTPCClustersArray *  GetCalcClustersArray(Bool_t newtree=kFALSE, const char* name=0)
  92 {
  93   //
  94   //construct AliTPCClusters Array object
  95   //
  96   AliTPCClustersArray * arr=0;
  97   //  if ( (gAlice!=0) && (gAlice->GetDetector("TPC")!=0) ) {
  98   //  arr = ((AliTPC*)gAlice->GetDetector("TPC"))->GetClustersArray();
  99   //  if (arr!=0) arr->Update();
 100   //}
 101   if (arr==0) {
 102     arr = new AliTPCClustersArray;
 103     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
 104     if (file==0) file = new TFile("galice.root","update");
 105     arr->SetClusterType("AliDigitCluster");
 106     arr->Setup(gTPCParam);
 107     cout<<"Update status : "<<arr->Update()<<"\n";
 108     char treeName[100];
 109     if (name==0)
 110       sprintf(treeName,"TreeCExact_%s",gTPCParam->GetTitle());
 111     else
 112       sprintf(treeName,"TreeCCalc_%s",gTPCParam->GetTitle());
 113     if (newtree!=kTRUE){
 114       cout<<"Connect tree status : "<<arr->ConnectTree(treeName)<<"\n";
 115     }
 116     else {
 117       arr->MakeTree();
 118     }
 119   }
 120   gCalcClusters = arr;
 121   return arr;
 122 }
 123
 124
 125 AliTPCClustersArray *  GetDifClustersArray(Bool_t newtree=kFALSE, const char* name=0)
 126 {
 127   //
 128   //construct AliTPCClusters Array object
 129   //
 130   AliTPCClustersArray * arr=0;
 131   //  if ( (gAlice!=0) && (gAlice->GetDetector("TPC")!=0) ) {
 132   //  arr = ((AliTPC*)gAlice->GetDetector("TPC"))->GetClustersArray();
 133   //  if (arr!=0) arr->Update();
 134   //}
 135   if (arr==0) {
 136     arr = new AliTPCClustersArray;
 137     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
 138     if (file==0) file = new TFile("galice.root","update");
 139     arr->SetClusterType("AliDifCluster");
 140     arr->Setup(gTPCParam);
 141     cout<<"Update status : "<<arr->Update()<<"\n";
 142     char treeName[100];
 143     if (name==0)
 144       sprintf(treeName,"TreeCExact_%s",gTPCParam->GetTitle());
 145     else
 146       sprintf(treeName,"TreeCDif_%s",gTPCParam->GetTitle());
 147     if (newtree!=kTRUE){
 148       cout<<"Connect tree status : "<<arr->ConnectTree(treeName)<<"\n";
 149     }
 150     else {
 151       arr->MakeTree();
 152     }
 153   }
 154   gDifClusters  = arr;
 155   return arr;
 156 }
 157
 158
 159
 160 TClonesArray * CompareClusters(TClonesArray *calcclusters,TClonesArray *exactclusters,
 161                                 TClonesArray *diffclusters, Float_t shx, Float_t shy)
 162 {
 163   if (calcclusters==0) return 0;
 164   if (exactclusters==0) return 0;
 165   if (diffclusters==0) return 0;
 166   TClonesArray & diff=*diffclusters;
 167
 168
 169   Int_t nclusters2 = calcclusters->GetEntriesFast();
 170   Int_t nclusters1 = exactclusters->GetEntriesFast();
 171
 172   for(Int_t i=0; i<nclusters1; i++){
 173     AliCluster * cl1 = (AliCluster*)exactclusters->UncheckedAt(i);
 174     AliDigitCluster *cl2;
 175     AliDifCluster diffc;
 176     Int_t index=-1;
 177     Float_t dx=10000;
 178     Float_t dy=10000;
 179     for (Int_t j=0; j<nclusters2;j++){
 180       cl2 = (AliDigitCluster*)calcclusters->UncheckedAt(j);
 181       if ( (cl2!=0)&& (cl1!=0)){
 182         Float_t ddx = (cl2->fX-shx)-cl1->fX;
 183         Float_t ddy = (cl2->fY-shy)-cl1->fY;
 184         if ((ddx*ddx+ddy*ddy)<(dx*dx+dy*dy)){
 185           dx=ddx;
 186           dy=ddy;
 187           index=j;
 188         }
 189       }
 190     }
 191
 192     cl2 = (AliDigitCluster*)calcclusters->UncheckedAt(index);
 193     if (cl2!=0){
 194       diffc.fDx      =dx;
 195       diffc.fDy      =dy;
 196
 197       diffc.fAngleX  = cl1->fSigmaX2;
 198       diffc.fAngleY  = cl1->fSigmaY2;
 199       diffc.fTracks[0] = cl2->fTracks[0];
 200       diffc.fTracks[1] = cl2->fTracks[1];
 201       diffc.fTracks[2] = cl2->fTracks[2];
 202       diffc.fGenerTrack= cl1->fTracks[0];
 203
 204       diffc.fX       = cl2->fX;
 205       diffc.fY       = cl2->fY;
 206       diffc.fNx      = cl2->fNx;
 207       diffc.fNy      = cl2->fNy;
 208       diffc.fQ       = cl2->fQ;
 209       diffc.fOrigQ   = cl1->fQ;
 210       diffc.fSigmaX2 = cl2->fSigmaX2;
 211       diffc.fSigmaY2 = cl2->fSigmaY2;
 212       diffc.fArea    = cl2->fArea;
 213       diffc.fMax     = cl2->fMax;
 214     }
 215     else cout<<"pici[ici/n";
 216
 217     new(diff[i]) AliDifCluster(diffc);
 218   }
 219   return diffclusters;
 220 }