MUON/MUONdigit.C

   1 void MUONdigit (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nCathode=1)
   2 {
   3 /////////////////////////////////////////////////////////////////////////
   4 //                                                                     //
   5 //   This macros converts the pad-hit information into digits          //
   6 //                                                                     //
   7 /////////////////////////////////////////////////////////////////////////
   8
   9 // Dynamically link some shared libs
  10
  11    if (gClassTable->GetID("AliRun") < 0) {
  12        gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so");
  13        gSystem->Load("$ALITOP/cern.so/lib/libPythia.so");
  14        gSystem->Load("$ROOTSYS/lib/libEG.so");
  15        gSystem->Load("$ROOTSYS/lib/libEGPythia.so");
  16        gSystem->Load("libGeant3Dummy.so");        //a dummy version of Geant3
  17        gSystem->Load("PHOS/libPHOSdummy.so");     //the standard Alice classes
  18        gSystem->Load("libgalice.so");             //the standard Alice classes
  19    }
  20
  21 // Connect the Root Galice file containing Geometry, Kine and Hits
  22
  23    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  24    if (file) file->Close();
  25    file = new TFile("galice.root","UPDATE");
  26    file->ls();
  27 //   file->Map();
  28
  29    printf ("\n File loaded !");
  30
  31 // Get AliRun object from file or create it if not on file
  32
  33    if (!gAlice) {
  34        gAlice = (AliRun*)file->Get("gAlice");
  35        if (gAlice) printf("\n AliRun object found on file\n");
  36        if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  37    }
  38    printf ("\n gAlice created !\n");
  39
  40    AliMUON *MUON  = gAlice->GetDetector("MUON");
  41
  42    AliMUONchamber*  iChamber;
  43    AliMUONsegmentation*  segmentation;
  44
  45    Int_t Npx[10];
  46    Int_t Npy[10];
  47
  48    Int_t nxmax=1026;
  49    Int_t nymax=1026;
  50    AliMUONlist *elem[10][1026][1026];
  51    TObjArray *obj=new TObjArray;
  52
  53    Int_t trk[500];
  54    Int_t chtrk[500];
  55
  56
  57    Int_t digits[3];
  58    TVector *trinfo;
  59    trinfo=new TVector(2);
  60
  61 //
  62 //   Loop over events
  63 //
  64
  65    Int_t Nh=0;
  66    Int_t Nh1=0;
  67    for (int nev=0; nev<= evNumber2; nev++) {
  68        Int_t nparticles = gAlice->GetEvent(nev);
  69        cout << "nev         " <<nev<<endl;
  70        cout << "nparticles  " <<nparticles<<endl;
  71        if (nev < evNumber1) continue;
  72        if (nparticles <= 0) return;
  73
  74        TTree *TH = gAlice->TreeH();
  75        Int_t ntracks = TH->GetEntries();
  76        Int_t Nc=0;
  77
  78        Int_t counter=0;
  79        //
  80        // loop over cathodes
  81        for (int icat=0; icat<nCathode; icat++) {
  82            //
  83            // initialize the array of pointers
  84            printf("\n Initialize array of pointers \n");
  85            for (int i=0; i<10; i++) {
  86                for (int j=0; j<nymax; j++) {
  87                    for (int k=0; k<nxmax; k++) {
  88                        elem[i][j][k]=0;
  89                    }
  90                }
  91            }
  92            printf("\n Finished !\n");
  93
  94 //
  95 //   Loop over tracks
  96 //
  97            printf("\n ntracks: %d",ntracks);
  98
  99            for (Int_t track=0; track<ntracks;track++) {
 100                gAlice->ResetHits();
 101                Int_t nbytes += TH->GetEvent(track);
 102                if (MUON)  {
 103                    for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1);
 104                        mHit;
 105                        mHit=(AliMUONhit*)MUON->NextHit())
 106                    {
 107                        Int_t   nch   = mHit->fChamber;  // chamber number
 108                        Float_t x     = mHit->fX;        // x-pos of hit
 109                        Float_t y     = mHit->fY;        // y-pos
 110
 111
 112                        iChamber = &(MUON->Chamber(nch-1));
 113                        response=iChamber->GetResponseModel();
 114
 115                        Int_t nsec=iChamber->Nsec();
 116                        Int_t rmin = (Int_t)iChamber->frMin;
 117                        Int_t rmax = (Int_t)iChamber->frMax;
 118                        if (nch > 10) continue;
 119
 120 //
 121 //
 122                        for (AliMUONcluster* mPad=
 123                                 (AliMUONcluster*)MUON->FirstPad(mHit);
 124                             mPad;
 125                             mPad=(AliMUONcluster*)MUON->NextPad())
 126                        {
 127                            Int_t cathode  = mPad->fCathode;   // cathode number
 128                            Int_t nhit     = mPad->fHitNumber; // hit number
 129                            Int_t qtot     = mPad->fQ;         // charge
 130                            Int_t ipx      = mPad->fPadX;      // pad  X
 131                            Int_t ipy      = mPad->fPadY;      // pad  Y
 132                            Int_t iqpad    = mPad->fQpad;      // charge per pad
 133                            Int_t rsec     = mPad->fRSec;      // sector#
 134 //
 135 //
 136                            if (cathode != (icat+1)) continue;
 137                            segmentation=iChamber->
 138                                GetSegmentationModel(cathode);
 139                            Int_t Npx[nch-1]  = segmentation->Npx();
 140                            Int_t Npy[nch-1]  = segmentation->Npy();
 141
 142
 143                            Int_t npx=Npx[nch-1];
 144                            Int_t npy=Npy[nch-1];
 145                            Float_t thex, they;
 146                            segmentation->GetPadCxy(ipx,ipy,thex,they);
 147
 148                            Float_t rpad=TMath::Sqrt(thex*thex+they*they);
 149
 150
 151                         // check boundaries
 152                            if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
 153
 154                            // fill the info array
 155
 156                            trinfo(0)=(Float_t)track;
 157                            trinfo(1)=(Float_t)iqpad;
 158
 159                            digits[0]=ipx;
 160                            digits[1]=ipy;
 161                            digits[2]=iqpad;
 162
 163                            // build the list of fired pads and update the info
 164
 165                            AliMUONlist *pdigit=elem[nch-1][ipy+npy][ipx+npx];
 166
 167                            if (pdigit==0) {
 168                                obj->AddAtAndExpand(
 169                                    new AliMUONlist(rpad,digits),counter);
 170                                counter++;
 171                                Int_t last=obj->GetLast();
 172                                pdigit=(AliMUONlist*)obj->At(last);
 173                                elem[nch-1][ipy+npy][ipx+npx]=pdigit;
 174                                // list of tracks
 175                                TObjArray *trlist=
 176                                    (TObjArray*)pdigit->TrackList();
 177                                trlist->Add(trinfo);
 178                            } else {
 179                                // update charge
 180                                (*pdigit).fSignal+=iqpad;
 181                                // update list of tracks
 182                                TObjArray* trlist=(TObjArray*)pdigit
 183                                    ->TrackList();
 184                                Int_t nptracks=trlist->GetEntriesFast();
 185
 186                                for (Int_t tr=0;tr<nptracks;tr++) {
 187                                    TVector *ptrk=(TVector*)trlist->At(tr);
 188 //                                 TVector &vtrk = *ptrk;
 189                                    Int_t entries=ptrk->GetNrows();
 190                                    trk[tr]=Int_t(ptrk(0));
 191                                    chtrk[tr]=Int_t(ptrk(1));
 192                                    if (trk[tr]==track) {
 193                                        chtrk[tr]+=iqpad;
 194                                        trlist->RemoveAt(tr);
 195                                        trinfo(0)=trk[tr];
 196                                        trinfo(1)=chtrk[tr];
 197                                        trlist->AddAt(trinfo,tr);
 198                                    } else {
 199                                        trlist->Add(trinfo);
 200                                    }
 201                                } // end loop over list of tracks for one pad
 202                            } // end if pdigit
 203                        } //end loop over clust
 204                    } // hit loop
 205                } // if MUON
 206            } // track loop
 207
 208            Int_t tracks[10];
 209            Int_t charges[10];
 210            cout<<"Start filling digits now !"<<endl;
 211            const Int_t zero_supm = 6;
 212            const Int_t adc_satm = 1024;
 213            Int_t nd=0;
 214
 215
 216            // start filling the digits
 217            for (Int_t id=0; id<10; id++) {
 218                Int_t nd=0;
 219                for (Int_t iy=0;iy<nymax;iy++) {
 220                    for (Int_t ix=0;ix<nxmax;ix++) {
 221                        AliMUONlist *address=elem[id][iy][ix];
 222                        if (address==0) continue;
 223                        // do zero-suppression and signal truncation
 224                        Int_t q=address->fSignal;
 225                        //
 226                        if ( q <= zero_supm ) continue;
 227                        if ( q > adc_satm)  q=adc_satm;
 228                        digits[0]=address->fPadX;
 229                        digits[1]=address->fPadY;
 230                        digits[2]=address->fSignal;
 231
 232                        TObjArray* trlist=(TObjArray*)address->TrackList();
 233                        Int_t nptracks=trlist->GetEntriesFast();
 234                        // this was changed to accomodate the real number of tracks
 235                        if (nptracks > 10) {
 236                            cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
 237                            nptracks=10;
 238                        }
 239                        for (Int_t tr=0;tr<nptracks;tr++) {
 240                            TVector *pp=(TVector*)trlist->At(tr);
 241                            Int_t entries=pp->GetNrows();
 242                            TVector &v2 = *pp;
 243                            tracks[tr]=(Int_t)v2(0);
 244                            charges[tr]=(Int_t)v2(1);
 245                        } //end loop over list of tracks for one pad
 246
 247                        // fill digits
 248                        MUON->AddDigits(id,tracks,charges,digits);
 249                        nd++;
 250                    } // end loop over ix
 251                } // end loop over iy
 252                cout<<"I'm out of the loops over pads"<<endl;
 253                cout<<" id  nd  "<<id<<" "<<nd<<endl;
 254            } //loop over chambers
 255            cout<<"I'm out of the loops for digitisation"<<endl;
 256
 257            gAlice->TreeD()->Fill();
 258            TTree *TD=gAlice->TreeD();
 259            int ndig=TD->GetEntries();
 260            cout<<"number of digits  "<<ndig<<endl;
 261
 262            TClonesArray *fDch;
 263            for (int i=0;i<10;i++) {
 264                fDch= MUON->DigitsAddress(i);
 265                int ndig=fDch->GetEntriesFast();
 266                printf (" i, ndig %d %d \n",i,ndig);
 267            }
 268
 269
 270            MUON->ResetDigits();
 271            obj->Clear();
 272        } //end loop over cathodes
 273
 274        char hname[30];
 275        sprintf(hname,"TreeD%d",nev);
 276        gAlice->TreeD()->Write(hname);
 277
 278        file->ls();
 279 //       file->Map();
 280    } // event loop
 281    obj->Delete();
 282    delete obj;
 283
 284    file->Close();
 285    cout<<"END  digitisation "<<endl;
 286
 287 }