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