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