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