ITS/AliITSFindClustersV2.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 /* $Id$ */
  17
  18 #include <TROOT.h>
  19 #include <TFile.h>
  20 #include <TTree.h>
  21 #include <TBranch.h>
  22 #include <TClonesArray.h>
  23 #include <TParticle.h>
  24
  25 #include "AliRun.h"
  26 #include "AliHeader.h"
  27
  28 #include "AliITS.h"
  29 #include "AliITSRecPoint.h"
  30 #include "AliITSFindClustersV2.h"
  31 #include "AliITSclusterV2.h"
  32 #include "AliITSgeom.h"
  33
  34 ClassImp(AliITSFindClustersV2)
  35
  36 //______________________________________________________________________
  37 AliITSFindClustersV2::AliITSFindClustersV2(){
  38     // Default constructor.
  39     // Inputs:
  40     //  none.
  41     // Outputs:
  42     //   none.
  43     // Return:
  44     //    A zero-ed constructed AliITSFindClustersV2 class.
  45
  46     fAr          = 0;
  47     fDeletfAr     = kFALSE; // fAr=0 dont't delete it.
  48     fGeom        = 0;
  49     fInFileName  = 0;
  50     fOutFileName = 0;
  51     fIn          = 0;
  52     fOut         = 0;
  53     fSlowFast    = kFALSE;  // slow simulation
  54     fInit        = kFALSE;  // Init failed
  55 }
  56 //______________________________________________________________________
  57 AliITSFindClustersV2::AliITSFindClustersV2(const TString infile,
  58                                            const TString outfile){
  59     // Standard constructor.
  60     // Inputs:
  61     //  const TString infile   Input file name where the RecPoints are
  62     //                         to be read from.
  63     //  const TString outfile  Output file where V2 tracking clulsters
  64     //                         are to be written. if =="" writen to the
  65     //                         same file as input.
  66     // Outputs:
  67     //   none.
  68     // Return:
  69     //    A standardly constructed AliITSFindClustersV2 class.
  70
  71     fAr          = 0;
  72     fDeletfAr    = kFALSE; // fAr=0 dont't delete it.
  73     fGeom        = 0;
  74     fInFileName  = 0;
  75     fOutFileName = 0;
  76     fIn          = 0;
  77     fOut         = 0;
  78     fSlowFast    = kFALSE;  // slow simulation
  79     fInit        = kFALSE;  // Init failed
  80
  81     fInFileName = new TString(infile);
  82     if(outfile.CompareTo("")!=0){
  83         fOutFileName = new TString(outfile);
  84     } // end if outfile.CompareeTo("")!=0
  85
  86     if(fOutFileName!=0){
  87         fIn = new TFile(fInFileName->Data(),"READ");
  88         fOut = new TFile(fOutFileName->Data(),"UPDATE");
  89     }else{ // open fIn file for update
  90         fIn = new TFile(fInFileName->Data(),"UPDATE");
  91     } // end if
  92
  93     fAr  = (AliRun*) fIn->Get("gAlice");
  94     if(!fAr){
  95         Warning("AliITSFindClusterV2",
  96                 "Can't fine gAlice in file %s. Aborting.",fIn->GetName());
  97         return;
  98     } // end if !fAr
  99     fDeletfAr = kTRUE; // Since gAlice was read in, delete it.
 100
 101     AliITS *its = (AliITS*) fAr->GetModule("ITS");
 102     if(!its){
 103         Warning("AliITSFindClusterV2",
 104                 "Can't fine the ITS in gAlice. Aborting.");
 105         return;
 106     } // end if !its
 107     fGeom = its->GetITSgeom();
 108     if(!fGeom){
 109         Warning("AliITSFindClusterV2",
 110                 "Can't fine the ITS geometry in gAlice. Aborting.");
 111         return;
 112     } // end if !fGeom
 113
 114     if(fOut) fOut->cd();
 115     fInit = kTRUE;
 116 }
 117 //______________________________________________________________________
 118 AliITSFindClustersV2::AliITSFindClustersV2(TFile *in,
 119                                            TFile *out){
 120     // Standard constructor.
 121     // Inputs:
 122     //  const TFile *in   Input file where the RecPoints are
 123     //                         to be read from.
 124     //  const Tfile *out  Output file where V2 tracking clulsters
 125     //                         are to be written. if ==0 writen to the
 126     //                         same file as input.
 127     // Outputs:
 128     //   none.
 129     // Return:
 130     //    A standardly constructed AliITSFindClustersV2 class.
 131
 132     fAr          = 0;
 133     fDeletfAr    = kFALSE; // fAr=0 dont't delete it.
 134     fGeom        = 0;
 135     fInFileName  = 0;
 136     fOutFileName = 0;
 137     fIn          = 0;
 138     fOut         = 0;
 139     fSlowFast    = kFALSE;  // slow simulation
 140     fInit        = kFALSE;  // Init failed
 141
 142     fIn  = in;
 143     fOut = out;
 144     fAr  = (AliRun*) fIn->Get("gAlice");
 145     if(!fAr){
 146         Warning("AliITSFindClusterV2",
 147                 "Can't fine gAlice in file %s. Aborting.",fIn->GetName());
 148         return;
 149     } // end if !fAr
 150     fDeletfAr = kTRUE; // Since gAlice was read in, delete it.
 151     AliITS *its = (AliITS*) fAr->GetModule("ITS");
 152     if(!its){
 153         Warning("AliITSFindClusterV2",
 154                 "Can't fine the ITS in gAlice. Aborting.");
 155         return;
 156     } // end if !its
 157     fGeom = its->GetITSgeom();
 158     if(!fGeom){
 159         Warning("AliITSFindClusterV2",
 160                 "Can't fine the ITS geometry in gAlice. Aborting.");
 161         return;
 162     } // end if !fGeom
 163
 164     if(fOut) fOut->cd();
 165     fInit = kTRUE;
 166 }
 167 //______________________________________________________________________
 168 AliITSFindClustersV2::AliITSFindClustersV2(AliRun *ar,
 169                                            const TString outfile){
 170     // Standard constructor.
 171     // Inputs:
 172     //  const TString infile   Input file name where the RecPoints are
 173     //                         to be read from.
 174     //  const TString outfile  Output file where V2 tracking clulsters
 175     //                         are to be written. if =="" writen to the
 176     //                         same file as input.
 177     // Outputs:
 178     //   none.
 179     // Return:
 180     //    A standardly constructed AliITSFindClustersV2 class.
 181
 182     fAr          = 0;
 183     fDeletfAr    = kFALSE; // fAr=0 dont't delete it.
 184     fGeom        = 0;
 185     fInFileName  = 0;
 186     fOutFileName = 0;
 187     fIn          = 0;
 188     fOut         = 0;
 189     fSlowFast    = kFALSE;  // slow simulation
 190     fInit        = kFALSE;  // Init failed
 191
 192     if(outfile.CompareTo("")!=0){
 193         fOutFileName = new TString(outfile);
 194     } // end if outfile.CompareeTo("")!=0
 195
 196     if(fOutFileName!=0){
 197         fOut = new TFile(fOutFileName->Data(),"UPDATE");
 198     } // end if
 199
 200     fAr  = ar;
 201     if(!fAr){
 202         Warning("AliITSFindClusterV2",
 203                 "ar==0. Aborting.");
 204         return;
 205     } // end if !fAr
 206     AliITS *its = (AliITS*) fAr->GetModule("ITS");
 207     if(!its){
 208         Warning("AliITSFindClusterV2",
 209                 "Can't fine the ITS in gAlice. Aborting.");
 210         return;
 211     } // end if !its
 212     fGeom = its->GetITSgeom();
 213     if(!fGeom){
 214         Warning("AliITSFindClusterV2",
 215                 "Can't fine the ITS geometry in gAlice. Aborting.");
 216         return;
 217     } // end if !fGeom
 218
 219     if(fOut) fOut->cd();
 220     fInit = kTRUE;
 221 }
 222 //______________________________________________________________________
 223 AliITSFindClustersV2::~AliITSFindClustersV2(){
 224     // Default constructor.
 225     // Inputs:
 226     //  none.
 227     // Outputs:
 228     //   none.
 229     // Return:
 230     //    A destroyed AliITSFindclustersV2 class.
 231
 232     fGeom = 0; // Deleted when AliRun/ITS is deleted.
 233     if(fInFileName!=0){ // input file opened by AliITSFindClustersV2
 234         if(fIn!=0){
 235             if(fIn->IsOpen()) fIn->Close();
 236             delete fIn;
 237             fIn = 0;
 238         } // end if
 239         delete fInFileName;
 240         fInFileName = 0;
 241     } // end if
 242
 243     if(fOutFileName!=0){ // input file opened by AliITSFindClustersV2
 244         if(fOut!=0){
 245             if(fOut->IsOpen()) fOut->Close();
 246             delete fOut;
 247             fOut = 0;
 248         } // end if
 249         delete fOutFileName;
 250         fOutFileName = 0;
 251     } // end if
 252     if(fDeletfAr && !fAr){
 253         cout << "deleting fAr."<< endl;
 254         delete fAr;
 255         fAr = 0;
 256         cout << "fAr deleted OK."<< endl;
 257     } // end if fDeletfAr
 258 }
 259 //______________________________________________________________________
 260 void AliITSFindClustersV2::Exec(const Option_t *opt){
 261     // Main FindclustersV2 function.
 262     // Inputs:
 263     //      Option_t * opt   list of subdetector to digitize. =0 all.
 264     // Outputs:
 265     //      none.
 266     // Return:
 267     //      none.
 268     Char_t name[50];
 269
 270     if(!fInit){
 271         Warning("Exec","(opt=%s) Initilization not succesfull. Aborting.\n",opt);
 272         return;
 273     } // end if !fInit
 274
 275     fGeom->Write();
 276
 277     fAr->GetEvent(0);
 278     TTree *pTree = fAr->TreeR();
 279     if(!pTree){
 280         Warning("Exec","Error getting TreeR. TreeR=%p",pTree);
 281         return;
 282     } // end if !pTree
 283     TBranch *branch = 0;
 284     if(fSlowFast) sprintf(name,"ITSRecPointsF");
 285     else sprintf(name,"ITSRecPoints");
 286     branch = pTree->GetBranch(name);
 287     if(!branch){
 288         Warning("Exec","Can't find branch %s in TreeR fSlowFast=%d",
 289                 name,fSlowFast);
 290         return;
 291     } // end if !branch
 292     TClonesArray *points = new TClonesArray("AliITSRecPoint",10000);
 293     branch->SetAddress(&points);
 294     Int_t nentr = (Int_t) branch->GetEntries();
 295
 296     if(fOut) fOut->cd();
 297     TClonesArray *cluster = new TClonesArray("AliITSclusterV2",10000);
 298     sprintf(name,"TreeC_ITS_%d",fAr->GetHeader()->GetEvent());
 299     TTree *cTree = new TTree(name,"ITS clusters");
 300     cTree->Branch("Clusters",&cluster);
 301     TClonesArray &cl = *cluster;
 302
 303     Float_t lp[5];
 304     Int_t lab[6];
 305     Int_t i,j,lay,lad,det,nclusters=0,ncl;
 306     Float_t kmip,x,y,zshift,yshift;
 307     Double_t rot[9];
 308     AliITSRecPoint *p;
 309     TParticle *part;
 310
 311     for(i=0;i<nentr;i++){
 312         points->Clear();
 313         branch->GetEvent(i);
 314         ncl = points->GetEntriesFast();
 315         if(ncl<=0) {
 316             cTree->Fill();
 317             continue;
 318         } // end if ncl<=0
 319         nclusters += ncl;
 320         fGeom->GetModuleId(i,lay,lad,det);
 321         fGeom->GetTrans(i,x,y,zshift);
 322         fGeom->GetRotMatrix(i,rot);
 323         yshift = x*rot[0] + y*rot[1];
 324         kmip = 1.0; // fGeom->GetModuletype(i)==kSPD
 325         if(fGeom->GetModuleType(i)==kSDD) kmip = 280.0;
 326         if(fGeom->GetModuleType(i)==kSSD) kmip = 38.0;
 327         for(j=0;j<ncl;j++){
 328             p = (AliITSRecPoint*)(points->UncheckedAt(j));
 329             lp[0] = - (p->GetX() + yshift);
 330             if(lay==1) lp[0] = -lp[0];
 331             lp[1] = p->GetZ() + zshift;
 332             lp[2] = p->GetSigmaX2();
 333             lp[3] = p->GetSigmaZ2();
 334             lp[4] = p->GetQ()/kmip;
 335             lab[0] = p->GetLabel(0);
 336             lab[1] = p->GetLabel(1);
 337             lab[2] = p->GetLabel(2);
 338             lab[3] = i;
 339             lad = lab[0];
 340             if(lad>=0){
 341                 part = (TParticle*) fAr->Particle(lad);
 342                 lad = -3;
 343                 while(part->P() < 0.005){
 344                     if(part->GetFirstMother()<0){
 345                         Warning("Exec","Primary momentum: %e",part->P());
 346                         break;
 347                     } // end if part->GetFirstMother()<0
 348                     lad = part->GetFirstMother();
 349                     part = (TParticle*) fAr->Particle(lad);
 350                 } // end while part->P() < 0.005
 351                 if(lab[1]<0) lab[1] = lad;
 352                 else if(lab[2]<0) lab[2] = lad;
 353                 else Warning("Exec","No empty lables!");
 354             } // end if lab>=0
 355             new(cl[j]) AliITSclusterV2(lab,lp);
 356         } // end for j
 357         cTree->Fill();
 358         cluster->Delete();
 359         points->Delete();
 360     } // end for i
 361     cTree->Write();
 362
 363     delete cTree;
 364     delete cluster;
 365     delete points;
 366 }